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Abstract 


\ 

Modern  optimal  control  methods  are  applied  to  a  lumped  mass 
model  of  a  tetrahedron.  The  four  unit  masses  of  this  model  are 
interconnected  by  isotropic  massless  rods  which  are  capable  of 
axial  deformation  only  (no  bending) .  NASTRAN  is  employed  in 
generating  a  normal  modes  approximation,  while  providing  the  mode 
shapes  and  frequencies  for  the  resultant  twelve  modes.  System 
control  is  achieved  via  collocated  sensor/actuator  pairs  at  three 
of  the  four  masses.  Pointing  accuracy  at  the  fourth  mass  is  used 
as  a  figure  of  merit  in  determining  the  effectiveness  of  the  con¬ 
troller.  A  prescribed  line  of  sight  response  is  established  as  a 
goal  for  successful  control. 

The  controller  is  developed  using  linear  optimal  techniques 
which  produce  feedback  gains  proportional  to  the  state.  The 
state  is  represented  as  modal  amplitudes  and  velocities  as  deter¬ 
mined  by  the  sensors.  The  four  higher  frequency  modes  are  trun¬ 
cated  to  signify  a  simplifying  order  reduction  step.  State  esti¬ 
mation  is  incorporated  due  to  the  non-availability  of  modal  ampli¬ 
tudes  and  velocities.  The  feedback  gains  are  established  via 
steady  state  optimal  regulator  theory.  Control  is  applied  with 
point  force  actuators.  System  response  is  examined  in  light  of 
the  effects  of  observation  spillover  and  control  spillover  onto 
a  specified  number  of  suppressed  modes.  v  A  comparison  is  ob¬ 
tained  by  complete  elimination  of  the  spillover  effect.  Using 
singular  value  decomposition,  the  spillover  is  first  eliminated 
through  judicious  reorientation  of  one  sensor/actuator  pair.  An 

attempt  to  control  two  modes  and  suppress  six  demonstrates  the 

vi 


advantages  of  spillover  elimination,  but  fails  to  satisfy  the 
specified  response  criteria.  /.  — 

Sensors  are  added  to  the  model  at  the  fourth  mass  and  ob¬ 
servation  spillover  is  again  eliminated.  Line  of  sight  response 
was  improved  over  the  case  without  sensor  additions,  but  was 
still  inadequate.  The  truncated  modes  were  added  to  the  system 
with  little  degradation,  verifying  the  acceptability  of  this 
truncation. 
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MODERN  OPTIMAL  CONTROL  METHODS 
APPLIED  IN  ACTIVE  CONTROL  OF  A 
TETRAHEDRON 

Introduction 

The  potential  for  larger  and  more  complex  space  structures 
has  grown  concurrent  with  the  approach  of  an  active,  operable 
space  transportation  program.  Present  system  concepts  involve 
the  deployment  of  earth  resource  satellites  and  micro-wave  power 
relay  systems  with  dimensions  extending  to  hundreds  and  eventually 
thousands  of  meters  in  diameter.  A  key  design  criteria  for  these 
immense,  mechanically  flexible  systems  is  the  requirement  to  de¬ 
velop  advanced  methods  for  control.  More  precisely,  a  principle 
issue  in  the  control  of  a  system  with  an  infinite  number  of  vibra¬ 
tional  modes  is  the  generation  of  a  method  for  stabilizing  these 
huge  structures  with  dimensionally  realistic  controllers.  This 
requirement  is  basically  a  function  of  on-board  computer,  sensor, 
and  force  actuator  limitations,  along  with  incumbent  modelling 
inaccuracies.  Of  the  numerous  methods  now  being  examined  as 
potential  solutions  to  this  control  problem,  modern  state  space 
control  theory  has  received  general  acceptance  as  the  most  via¬ 
ble  technique.  Applying  classical  control  methods  to  these  large 
structures  is  seen  as  computationally  improbable;  at  the  same 
time,  modern  state  space  theory,  incorporating  a  finite  element 
system  representation,  can  be  successfully  applied  to  a  very  wide 


class  of  flexible  structures.  This  theory  is  most  commonly 
applied  using  an  optimal  time- invariant  linear  regulator  as  a 
means  of  actively  controlling  vibration. 

Due  to  the  inherent  hardware  limitations  briefly  highlighted, 
active  control  must  be  restricted  to  a  relatively  small  number 
of  critical  modes.  Therefore,  in  a  necessary  truncation  step, 
some  higher  frequency  modes  remain  unmodelled.  Natural  damping 
in  the  system  is  assumed  to  preclude  the  possiblity  of  instabil¬ 
ity  resulting  from  these  modes.  Of  the  remaining  modes  (still  a 
potentially  large  number)  it  is  further  desirable  to  treat  only 
a  critical  few  (not  necessarily  those  with  the  lowest  frequency) , 
while  suppressing  the  rest. 

However,  the  sensor  outputs  are  contaminated  by  the  remaining 
"suppressed"  modes,  and  the  eventual  feedback  control  also  excites 
these  modes.  Balas  (Ref  1)  labels  these  effects  "observation 
spillover"  and  "control  spillover"  respectively.  He  shows  that 
either  or  both  of  these  effects  can  lead  to  overall  instabilities; 
the  suppressed  modes  must,  as  a  result,  be  a  design  consideration. 
Balas  describes  a  technique  with  which  to  develop  a  feedback  con¬ 
troller  using  state  variable  methods.  The  key  to  this  approach  is 
the  use  of  narrow  bandpass  filters  which  effectively  comb  out 
the  suppressed  mode  frequencies  to  eliminate  observation  spillover. 

Another  method  for  developing  an  appropriate  feedback  con¬ 
troller  was  first  presented  by  Sesak  (Ref  2) ,  and  later  expanded 
by  Coradetti  (Ref  3) .  This  approach  involves  the  use  of  a  so- 
called  "singular  perturbation"  technique  in  analyzing  and  elim¬ 
inating  the  spillover-generated  instabilities.  Coradetti  concludes 
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that  employing  this  "singular  perturbation"  method  in  a  limiting 
sense,  with  an  infinite  penalty  applied  against  any  spillover,  is 
equivalent  to  finding  a  transformation  matrix.  This  transforma¬ 
tion  matrix,  when  applied  to  feedback  gains,  effectively  elimi¬ 
nates  any  spillover  terms.  It  should  be  noted  that,  even  if 
spillover  does  not  render  the  system  unstable,  applying  the  trans¬ 
formation  method  may  still  improve  performance.  Additionally, 
while  no  method  for  actually  automating  optimal  sensor  and  actua¬ 
tor  placement  is  defined,  some  valuable  insight  into  the  nature 
of  this  task  is  precipitated.  This  is  accomplished  utilizing 
what  have  now  become  well  known  state  space  control  techniques 
in  conjunction  with  singular  value  decomposition  of  the  rectan¬ 
gular  matrices  of  modal  amplitudes  (Ref  4) . 

The  principle  function  of  this  thesis  is  to  provide  appli¬ 
cation  of  the  Coradetti  approach  to  a  three  dimensional,  lumped 
mass  model  of  a  tetrahedron.  A  line  of  sight  at  one  of  the  masses 
(simulating  pointing  accuracy)  will  be  used  as  a  figure  of  merit 
with  which  to  judge  the  general  effectiveness  of  this  method. 

This  thesis  will  serve  as  a  direct  extension  of  the  work  done  by 
Sanborn  (Ref  5) ,  in  which  the  stability  of  a  cantilever  beam  in 
bending  vibration  was  studied.  Specifically,  this  thesis  will 
examine  model  response  as  affected  by  the  number  and  orientation 
of  position  sensors  and  force  actuators.  The  elimination  of 
control  spillover  and  observation  spillover  will  be  obtained  using 
singular  perturbation  and  singular  value  decomposition  techniques. 

A  representation  of  a  tetrahedron  has  been  obtained  via 
the  normal  modes  approximation  package  found  in  the  NASTRAN  finite 
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element  computer  program.  The  natural  frequencies  and  eigen¬ 
values/eigenvectors  associated  with  each  mode  were  provided  by  a 
study  done  by  the  Charles  Stark  Draper  Laboratory.  For  appli¬ 
cation  of  the  control  method,  position  sensors  are  used  to  eval¬ 
uate  modal  amplitudes,  while  point  force  actuators  accomplish 
the  state  variable  feedback  control.  Singular  value  decomposi¬ 
tion  of  the  matrices  of  modal  amplitudes  at  sensor  locations 

/  4 

and  actuator  locations  is  used  to  produce  a  transformation  matrix 
by  which  spillover  terms  are  eliminated.  A  model  with  higher 
order  modes  truncated  (un-modelled)  is  used  to  design  the  control¬ 
ler.  The  effectiveness  of  this  controller  against  all  modes  is 
examined.  Finally,  a  study  of  improved  performance  with  added 
sensors  is  generated. 
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System  Model 


General  Configuration 

Of  the  many  design  criteria  which  must  be  considered  for 
the  large  flexible  spacecraft  currently  being  advanced,  pointing 
accuracy  looms  as  the  most  critical.  As  a  function  of  system 
size  and  operating  frequency,  pointing  accuracies  in  the  range 
of  one  tenth  the  half-power  signal  beam  width  will  be  required. 
The  ability  to  meet  these  stringent  requirements  becomes  a  direct 
function  of  the  isotropic  stiffness  of  the  system.  One  of  the 
space  erectable  or  assembly  concepts  that  has  the  promise  of 
supplying  this  needed  stiffness  in  larger  systems  is  the  geo¬ 
detic  truss  (Ref  6) .  Based  on  current  Space  Shuttle  cargo 
capacities,  whole  units  of  up  to  91.4  meters  can  be  packaged  for 
deployment.  For  very  large  systems,  these  units  are  assembled 
as  an  amalgamation  of  tetrahedrons —  the  basic  unit  of  geodetic 
truss.  By  changing  the  size  of  the  tetrahedrons,  a  large  array 
of  varying  stiffness  antenna  substructures  can  be  developed. 

For  this  reason,  a  tetrahedron  is  seen  as  an  important  model 
against  which  to  apply  proposed  control  techniques. 

The  finite  element,  lumped  mass  model  to  be  used  herein 
is  depicted  in  Fig  1  and  Fig  2.  This  model  is  seen  to  consist 
of  ten  nodes.  The  twelve  interconnecting  truss  members  are 
assumed  to  be  massless  and  are  capable  of  resisting  or  exerting 
axial  force  only  (no-bending).  Masses  are  of  one  unit  each, 
and  are  located  at  grid  points  one  through  four.  Each  mass 
is  capable  of  perturbation  with  three  degrees  of  freedom. 
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The  remaining  grid  points  (five  through  ten)  serve  to  esta¬ 
blish  a  fixed  line  of  sight  for  an  initial  set  of  six  col¬ 
located  sensor /actuator  pairs.  Node  coordinates  for  the  model 
are  listed  as  Table  I.  For  this  analysis,  position  sensors 
are  employed,  but  velocity  sensors  are  not.  The  effects  of 
this  are  detailed  in  the  linear  system  model  to  be  developed 
in  this  section. 

An  eigenvalue  analysis  of  this  nominal  model  has  been  pro¬ 
vided  via  the  NASTRAN  computer  program.  Key  results  of  this 
analysis  are  presented  as  Table  II.  The  eigenvectors  associated 
with  these  eigenvalues  can  be  found  in  Appendix  A.  Table  III 
provides  the  initial  conditions  required  for  a  time  history  exam¬ 
ination  of  system  stability.  This  stability  will  be  assessed 
using,  as  a  figure  of  merit,  the  pointing  accuracy  along  the  Z 
axis  at  node  1.  Since  any  perturbation  directly  along  the  Z  axis 
has  no  impact  on  pointing  accuracy,  the  line  of  sight  in  the 
X  and  Y  directions  only  will  be  examined. 


Table  I. 

Node  Coordinates 


Node 

X 

Y 

Z 

1 

0.0 

0.0 

10.165 

2 

-5.0 

-2.887 

2.0 

3 

5.0 

-2.887 

2.0 

4 

0.0 

5.7735 

2.0 

5 

-6.0 

-1.1547 

0.0 

6 

-4.0 

-4.6188 

0.0 

7 

4.0 

-4.6188 

0.0 

8 

6.0 

-1.1547 

0.0 

9 

2.0 

5.7735 

0.0 

10 

-2.0 

5.7735 

0.0 
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Table  II 


Key  Results  of  NASTRAN  Eigenvalue  Analysis 


Mode 

Generalized 

Mass 

Generalized 

Stiffness 

u)n(rad)/8ec 

ft  (rad 

1 

1.0E+00 

1.37E+00 

1.17E+00 

1.37E+00 

2 

1.0E+00 

2.15E+00 

1.47E+00 

2.15E+00 

3 

1.0E+00 

8.79E+00 

2.96E+00 

8.79E+00 

4 

1.0E+00 

1.26E+01 

3.56E+00 

1.26E+01 

5 

1.0E+00 

1.48E+01 

3.85E+00 

1.48E+01 

6 

1.0E+00 

2.65E+01 

5.15E+00 

2.65E+01 

7 

1.0E+00 

3.22E+01 

5.67E+00 

3.22E+01 

8 

1.0E+00 

3.26E+01 

5.71E+00 

3.26E+01 

9 

1.0E+00 

7.99E+01 

8.93E+00 

7 . 99E+01 

10 

1.0E+00 

1.06E+02 

1.03E+01 

1.06E+02 

11 

1.0E+00 

1.19E+02 

1.09E+01 

1. 19E+02 

12 

1.0E+00 

1.95E+02 

1.40E+01 

1.95E+02 

Table  III 

Initial 

Conditions  for  Time  History 

Response 

Mode 

Displacement  (n) 

Velocity  (n) 

1 

-.001 

-.003 

2 

.006 

.01 

3 

.001 

.03 

4 

-.009 

-.02 

5 

.008 

.02 

6 

-.001 

-.02 

7 

-.002 

-.003 

8 

.002 

.004 

9 

.0 

.0 

10 

.0 

.0 

11 

.0 

.0 

12 

.0 

.0 
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Equations  Of  Motion 

Since  there  are  no  exact  equations  of  motion  for  a  con- 
tinous  model  of  a  tetrahedron,  we  are  restricted  to  the  discre¬ 
tized  representation  provided  by  the  finite  element  routines. 

The  output  function  or  motion  of  the  model  can  be  expressed  as: 

n 

Y(x.,t)  =  Z  <J> .  (x  . )  U .  (t)  (1) 

J  i=l  J  i 

where  the  <Jk  (x^)  terms  are  the  mode  shapes,  and  the  IL  (t)  terms 
are  the  mode  amplitudes,  with  n  being  the  number  of  modes  ex¬ 
hibited  by  the  model.  For  an  exact  solution  to  a  continuous 
system,  the  number  of  lumped  masses  and  the  number  of  modes  (n) 
would  have  to  reach  infinity.  Practically  speaking,  the  total 
system  displacement  Y(Xj,t),  can  be  reasonably  represented  by  a 
truncation  of  Eq  (1)  such  that  n  is  considerably  less  than  infin¬ 
ity.  This  truncation  will,  of  course,  lead  to  model  reduction 
errors;  but,  up  to  a  certain  point  these  errors  are  relatively 
insignificant. 

NASTRAN  analyzes  the  model  in  Fig  1.  and  generates  both 
the  normal  mode  shapes  and  the  corresponding  natural  frequencies 
(u>n)  .  Since  this  is  a  lumped  mass  model  consisting  of  four  masses, 
with  each  mass  having  three  degrees  of  freedom,  there  are  a  total 
of  twelve  normal  modes. 

Linear  System  Model 

As  stated,  the  number  of  modes  (n)  for  a  complex  model  may 
be  very  large.  The  practical  limitations  for  an  on-board  com¬ 
puter  and  the  associated  sensor  and  actuator  hardware  make  it 
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necessary  to  develop  a  controller  that  is  concerned  with  a 
minimum  number  of  modes,  while  still  satisfying  what  may  be  very 
stringent  requirements  on  the  performance  (here,  line  of  sight 
accuracy) .  As  the  control  theory  outlined  in  this  paper  is  elab¬ 
orated,  a  possible  method  for  determining  which  modes  require 
control  will  be  discussed.  At  this  point,  assuming  this 
determination  is  possible,  the  system  output  of  Eq  (1)  can  be 
segregated  into  3  partitions;  the  controlled,  the  suppressed, 
and  the  unmodelled: 

=  Yc(xj,t>  +  Ys(Xj,t)  +  Y^IXjt)  (2) 

Yum(Xj,t)  is  that  portion  of  the  output  generated  through 
the  highest  frequency  modes.  These  modes  are  unmodelled,  with 
the  hope  that  the  bandwidths  of  the  sensors  and  actuators  employ¬ 
ed  will  be  less  that  the  natural  frequencies  of  the  modes.  Fur¬ 
thermore,  since  these  modes  have  such  high  frequencies,  they  may 
be  quite  difficult  to  excite.  Hence,  any  controller  designed 
for  this  system  can  ignore  these  modes.  These  modes  are  sub¬ 
sequently  called  the  residual  modes. 

Y  (x .  t)  is  that  portion  of  the  output  generated  by  modes 
s  j , 

of  less  high  frequency,  which,  none-the-less  have  a  minimal 
direct  impact  on  system  performance.  Due  to  their  indirect 
and  potentially  destabilizing  impact  (spillover) ,  they  must 
be  included  in  the  design  process.  These  modes  are  subsequently 
called  the  suppressed  modes. 

Yc(Xj,t)  is  that  portion  ot  the  output  which  we  must  directly 
control  to  insure  satisfactory  performance.  These  critical  modes 


will  subsequently  be  called  the  controlled  modes. 

Equation  (2)  can  now  be  written  in  segregated  form  as: 

c 


Y  (x. ,t)  =  Z  <J>. (x  . ) U .  (t) 
c  3  i=l  1  3  i 

(3) 

c+s 

Y  (x.  ,t)  =  E  4> .  (x  . )  U .  (t) 

s  3  i=c+l  1  3  x 

(4) 

n. 

(5) 

where  c  is  the  number  of  controlled  modes,  s  is  the  number  of 
suppressed  modes,  and  n  is  the  total  number  of  modes  in  the 
model.  Again,  for  this  system  model,  n  is  twelve.  For  the  pur¬ 
pose  of  future  analysis,  the  case  of  truncating  the  highest  fre¬ 
quency  modes  will  be  simulated  by  suggesting  that  the  last  four 
(highest  natural  frequency)  modes  generated  by  NASTRAN  fall  into 
this  category.  The  design  process  for  the  overall  controller 
will  be  based  on  knowledge  of  only  the  first  eight  modes.  The 
eventual  controller  will  be  applied  to  a  system  incorporating 
all  twelve  modes  in  an  attempt  to  verify  the  acceptability  of 
this  truncation.  The  modelling  can  thus  be  seen  as  a  process  of 
two  truncations  in  the  effort  to  reduce  control  hardware  and 
software  requirements.  First,  the  model  is  truncated  to  a  work¬ 
able  number  of  modes  by  designing  a  controller  that  is  blind  to 
the  higher  frequency  modes.  Second,  the  model  is  limited  to 
controlling  only  the  critical  modes  where  the  figure  of  merit 
is  concerned. 

NASTRAN  has  taken  the  prescribed  system  with  the  masses 
and  gridpoints  provided,  and  modelled  the  structure  with  a  set 


12 


of  second  order  differential  equations.  These  are  the  basic 


spring  mass  differential  equations  such  that  ri  +  =  f. 

The  associated  first  order  eigenproblem  is  solved  (Ref  7)  so 
as  to  provide  the  decoupled  normal  modes.  This  allows  assem¬ 
bling  a  state  space  representation  of  the  system: 

X(t)  =  AX(t)  +  Bu  (t)  (6) 


where 


X(n  x  1)  is  the  state  vector 

u(m  x  1)  is  the  control  input  vector 

A(n  x  n)  is  the  plant  matrix 

B(n  x  m)  is  the  input  matrix 

By  letting  the  state  X  be  the  partitioned  matrix  of  mode  ampli¬ 
tudes  (U^t))  and  their  rates  of  change  (U.(t))  the  state  variables 
become : 


X(t) 


— T 

It 

ui(t) 

tK(t) 

j 

i  1,2, ... ,n  (7) 


Further  separating  the  states  into  X  ,  formed  by  the 

c 

controlled  amplitudes  and  rates;  and  X  ,  formed  by  the  suppressed 

s 

amplitudes  and  rates  renders: 


Xc(t)  -  | 

u±  (t) 

crAct|  | 

I  T 

i  —  1 , 2 , . .  ,  c 

(8) 

v« ' 

Uj(t) 

Uj  (t) 

T 

j  =  c+1 , . .c+s 

(9) 

Substituting  these  states  into  Eq  (6) ,  the  system  is  now  modelled 
by: 


X  (t)  =  A  X  (t)  +  B  u  ( t )  (10) 

c  c  c  c 

X_(t)  =  A  X  (t)  +  B  u ( t)  (11) 

s  s  s  s 
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The  system  parameter  matrices  are  defined  as: 


(12) 


(13) 


(14) 


(15) 


The  uL  and  ok  terms  are  the  diagonal  elements  of  square 
matrices  which  represent  the  natural  frequencies  of  the  con¬ 
trolled  and  suppressed  modes  respectively,  while  the  and 
5  terms  represent  the  damping  ratios  for  those  modes  i  =  l,2,..c 

and  j  =  c+l,c+2, . . .c+s.  The  ?2  and  ft  terms  are  diagonal  matrices 

c  s 

of  these  same  natural  frequencies  squared  as  determined  by 
NASTRAN.  Therefore,  as  an  example,  if  two  modes  for  a  given 
system  were  to  be  controlled,  one  would  have: 
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0  0  10 

0  0  0  1 

A  =  ( 

c  -w2  0  -2K  03  0 

l  ii 

0  -to2  0  -25  u> 

2  2  2 

Furthermore,  the  B  and  B  matrices  aethe  control  input 
matrices,  and  are  those  matrices  whose  columns  are  the  mode 
shapes  (<p^  (x) ,  <p_.  (x) )  evaluated  at  each  actuator  location  such 
that: 


B  = 
-c 


(J>1(xl)  2) 

4>2(x1)  4>2 


•  •  4*2 
■  ' ^2 


B  = 


c(xl> 

♦c'x2>  •  • 

'*e(xa,J 

c+l(xl> 

♦c+l'V. 

.  .WXa> 

C+2(X1> 

• 

Wx2> 

^c+2  txa) 

• 

c+s|xl» 

*c*s,x2>- 

*  •U(xa) 

(17) 


(18) 


where  a  is  the  total  number  of  actuators  employed. 

Additionally,  state  space  methods  render  the  sensor  output 

as: 


Y(t)  =  CcXc(t)  +  CsXs(t) 


(19) 


with 


(20) 
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c 

-  s 


0 


(21) 


c  = 

a 


j 


where  C  and  C  are  matrices  whose  "rows"  are  the  mode  shapes  of 

W  S 

the  controlled  and  the  suppressed  modes  respectively  evaluated 
at  the  prescribed  sensor  locations  such  that: 

(22) 


(23) 

The  null  portion  of  the  C  and  C  matrices  represent 

~  c  ~  s 

the  velocities  at  the  prescribed  sensor  locations,  which  are 
zero  since  only  displacement  sensors  are  being  employed.  Again, 
b  is  the  total  number  of  sensors  used.  It  should  be  clear  that 
if  collocated  sensors  and  actuators  are  used,  with  a  =  b,  then 


C  = 
~c 


*l(xl>  *2<xl) 
^ 1 ^X2 ^  ^2  ^x2 ^ 


w  w 


c  = 


*c+2(V 


W 

VX2> 


W 


0+1 <xl' 

*C+2(X1K 

•  •  *c+s(xl) 

C+1<X2> 

4’o+2(X2l- 

•  •  (f)c+s(x2) 

Wxb> 


B  T  =  C  (24) 

~c  ~  c 

BcT  =  C  (25) 

~  s  ~s 


As  these  model  elements  are  created,  it  becomes  clear  that 
this  methodology  is  independent  of  structural  complexity,  except 
for  the  overall  matrix  dimensions.  Therefore,  the  applicability 
of  the  subsequent  analysis  can  be  seen  to  be  far  reaching. 


1 
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As  a  starting  point  toward  developing  a  state  variable 
feedback  controller.  Fig  3  below  represents  the  uncontrolled 
system  that  has  been  here-to-fore  described. 


Figure  3.  Simple  Open  Loop  Plant 

In  order  to  eventually  form  an  active  control,  u(t),  using 
state  variable  (modern  control)  feedback  techniques,  complete 
knowledge  of  the  actual  state  at  time  t  must  be  known.  However, 
the  only  measure  of  the  state  X  is  the  measurement  vector  Y  pro¬ 
vided  by  the  sensors.  To  take  those  observations  and  create  the 
corresponding  state,  it  will  be  necessary  to  develop  a  state  esti¬ 
mator  which  will  accept  those  sensor  observations  and  estimate 
X  as  X. 

Modal  Control 

As  Balas  explains,  the  state  estimator  used  in  developing 
active  feedback  control  can  either  be  a  Kalman  Filter  when  it 
is  found  that  the  signal-to-noise  ratios  are  relatively  small, 
or  a  Luenberger  observer,  or  a  least  squares  technique,  should 
the  signal  to  noise  ratio  be  high  enough  to  treat  the  system 
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as  deterministic.  Regardless  of  which  is  used,  the  estimator 
will  have  the  form; 


Vfc)  =  ¥n  +  V(t)  +  KN[Y(t)  “  Y(tn  (26) 

and  Xx7(0)  =  0  (27) 

N 

Y(t)  =  CNXN  (28) 

where  N  is  replaced  in  our  system  by  either  c  or  s\ 

Observation  of  Eq  (26)  shows  that  the  estimator  equation 
is  comprised  of  the  internal  model  of  the  state  as  in  Eqs  (10) 
and  (11)  ,  plus  a  correction  term  which  is  made  up  of  the  error 
between  the  measured  output  (Y(t))  and  the  computed  output 
(Y (t) ) .  Equation  (27)  establishes  an  initial  condition  for  the 
state  out  of  convenience.  The  error  in  this  state  estimation 
process  is  given  as: 

eN(t>  =  Vt)  -  XN(t)  (29) 

The  equations  for  this  estimator  error,  formed  by  combining 
Eqs  (26),  (27),  and  (28)  with  Eqs  (10),  (11),  and  (19)  becomes: 

iN(t)  =  (^  -  KNCK);s(t.|  +  KNCRXR(t)  (30) 

For  the  prescribed  system,  this  becomes: 

e  (t)  =  (A  -  K  C  |e  It)  +  K  C  x  (t)  (31) 

c  c  ccc  css 

Ignoring  the  suppressed  modes,  this  finally  becomes: 

ic(t)  -  <AC  -  KcCc)5c(t)  (32) 
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The  observer  gain  matrix,  K,  must  be  formulated  so  as  to  in¬ 
sure  that  the  estimator  error  defined  in  Eq  (32)  decays  expon¬ 
entially  at  a  rate  more  rapid  than  the  system  dynamics.  The 


decay  rate  is  determined  by  the  eigenvalues  of  (A  -KC  ) .  Since 

c  c 

the  eigenvalues  of  a  matrix  are  equal  to  the  eigenvalues  of 
the  transpose  of  that  matrix,  Eq  (32)  can  be  rewritten  as 
follows: 


w(t)  = 


A  ^(t)  - 
c 


ccT9(t) 


(33) 


g  (t)  =  KTw  ( t ) 


(34) 


The  observer  gain  matrix,  K,  can  now  be  calculated  via 
steady  state  optimal  regulator  theory.  This  is  equivalent  to 
minimizing  the  quadratic  regulator  performance  index  J,  where: 

CO  _T  _  — fP  _ 

J  =  Hf  o  (wQobw  +  g^R^gJdt  (35) 


The  known  optimal  solution  to  this  minimization  problem  is: 


K 


T 


CcP 


(36) 


where  P  is  the  solution  to  the  steady  state  algebraic  matrix 
Ricatti  Equation: 

PA  +  A  P  -  PC  TR  P  +  Q  .  =0  (37) 

c  c  c  ob  c  ob 

where  Q  .  is  an  n  x  n  positive  semidefinite  state 

weighting  matrix 

R  ,  is  an  m  x  m  positive  definite  control 
weighting  matrix 

By  treating  only  the  controlled  modes  in  the  generation 
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of  the  optimal  state  feedback  gain  matrix,  we  have  signifi¬ 
cantly  reduced  the  order  of  the  controller.  This  was  accom¬ 
plished,  as  previously  stated  to  avoid  practical  problems 
encountered  in  deriving  a  global  controller.  The  reduced  or¬ 
der  controller  will  subsequently  be  designed  to  control  a  sub¬ 
set  of  all  of  the  system  states,  while  simultaneously  avoiding 
any  excitation  of  the  remaining  states.  Coradetti  clarifies 
the  advantages  of  this  process  when  he  points  out  that  the 
computational  burden  of  solving  the  Ricatti  Equation  increases 
roughly  as  the  cube  of  the  order  of  the  equation.  There  may 
simply  not  be  sufficient  on-board  computer  memory  available. 

Also,  the  state  estimator  process  increases  with  system  order  a 
at  a  greater  than  linear  rate.  Finally,  with  non-interacting 
controller  there  will  be  greater  fault  tolerance  to  actuator 
failures . 

In  precisely  the  same  fashion,  the  control  feedback  gain 
matrix,  G,  can  be  formulated.  Now,  again  using  steady  state 
optimal  regulator  theory,  the  performance  index  to  be  minimized 
is: 

J  =  hf°°(X  TFX  +  fTRf)dr  (38) 

0  c  c 

where  F  is  an  n  x  n  positive  semidefinite  state 

weighting  matrix 

R  is  an  m  x  m  positive  definite  control 
weighting  matrix 

The  optimal  solution  to  this  minimization  problem  is 


where  s  is  the  solution  to  the  matrix  Ricatti  Equation: 


SA  +  A  TS  -  SB  R_1B  TS  +  F  =  0  (40) 

Implementing  the  results  of  Eqs  (36)  and  (39)  with  the 
system  Eqs  (10)  and  (11),  as  well  as  Eq  (32)  renders: 

» 

X  (t)  =  (A  +  B  G ) X  ( t )  +  B  Ge(t)  (41) 

v  c  c  c  c 

Xs(t)  =  AsXs(t)  +  BsGXc(t)  +  BgGe (t)  (42) 


By  taking  the  step  of  defining  a  system  state  vector  incor¬ 
porating  the  controlled  states,  the  suppressed  states,  and 
the  estimator  error,  such  that: 


Z(t) 


4 


—  T 

Xc  (t) 


— T 

eUt) 


—  T 

xs  (t) 


T 


(43) 


a  closed  loop  system  model,  containing  the  effects  of  the 
suppressed  and  controlled  modes,  and  utilizing  state  variable 
feedback  as  the  control  mechanism  can  be  presented  as: 


Z(t)  = 


A  +B  G 
c  c 

bcg  : 

0 

0 

u 

1 

< 

KC 

c  c 

s 

B  G 

B  G 

A_ 

S 

s 

S 

(44) 


Recalling  that  the  observation  and  control  feedback  gain 

matrices  (K  and  G)  were  designed  to  operate  on  the  controlled 

modes,  the  terms  KC  and  B  G  create  potential  problems.  These 

s  s 

in  effect,  are  known  as  observation  spillover  and  control  spill¬ 
over,  respectively.  Although  all  of  the  diagonal  matrices  of 


Eq  (44)  are  designed  to  have  purely  negative  real  parts  for 

all  eigenvalues,  it  is  obvious  that  the  KC  and  B  G  terms  can 

s  s 

cause  overall  system  instabilities. 


Block  Diagram  Representation  for  the  Linear  Model 


In  a  parallel  section  of  his  paper,  Sanborn  generates  the 
block  diagrams  representing  this  new  system  model  in  two  seper- 
ate  forms.  Since  the  equations  now  governing  the  model  are: 


X  =  AX  +  Bu 


State  Equation 


(45) 


Y  =  CX 


Output  Equation  (46) 


u  =  GXc  Control  Equation  (47) 

A  A.  A 

&  =  A  X  +Bu+K(Y-Y)  Estimator  Equation  (48) 

C  CO  c 

The  system  can  be  presented  as  Fig  4.  This  diagram  can 
be  manipulated  per  Johnson  (Ref  8)  to  generate  a  modified  block 
diagram  form  as  shown  in  Fig  5. 

From  Fig  5  the  closed  loop  transfer  function  for  the  con¬ 
troller  is  seen  to  be: 

J(S)  =  K ( SI  -  A  -B  G  +  KC  ) _1G  (49) 

y  c  c  c 

From  this  transfer  function,  we  know  that  if  any  of  the 

eigenvalues  of  (A  +  B  G  -  KC  )  are  positive,  then  the  controller 

c  c  c 

is  unstable.  Since  the  techniques  for  generating  both  the  ob¬ 
servation  and  control  gain  matrices  were  employed  independently, 
the  possibility  that  an  unstable  controller  is  formed  exists. 
Although  the  controller,  when  coupled  with  the  plant,  would 
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Plant 


Figure  4.  System  in  Block  Diagram  Form 

Plant 


Controller 


Figure  5.  System  in  Modified  Block  Diagram  Form 
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9 


produce  a  stable  system,  the  potentially  disastrous  effects  of 
an  intermittent  decoupling  must  be  emphasized.  An  examination 
of  the  eigenvalues  of  (Ac  +  BcG  -  KCc)  will,  therefore,  be 
included  in  the  analysis. 
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Transformation  Matrix  Control 


It  has  been  shown  in  what  has  preceded  that,  due  to  obser¬ 
vation  spillover  and  control  spillover,  the  system  represented 
by  Eq  (44)  could  be  made  unstable.  In  an  attempt  to  alleviate 
this  problem  we  will  employ  a  control  technique  which  attempts 
to  eliminate  spillover.  This  suggests  driving  the  off-diagonal 
matrices  of  Eq  (44)  to  zero,  while  retaining  active  feedback 
control  of  the  overall  system.  An  examination  of  the  system 

equation  leads  one  to  realize  that,  if  either  B  G  or  KC  are 

s  s 

zero,  the  system  eigenvalues  revert  to  the  eigenvalues  falling 

on  the  diagonal.  The  nature  of  these  diagonal  matrices  is 

such  that  negative  eigenvalues  (and,  hence,  system  stability) 

are  guaranteed.  Obviously,  one  solution  to  B  G  =  0  is  G  =  0. 

s 

However,  this  solution  also  renders  B^G  =  0,  and  control  is  for¬ 
gone.  That  being  the  case,  the  transformation  method  is  dir¬ 
ected  at  constraining  the  feedback  gain  matrices  such  that: 

B  G  =  0 
s 

KC  =  0 
s 

while,  at  the  same  time: 


B  G  #  0 

(52) 

c 

KC.  *  0 

(53) 

To  develop  this  method,  we  will  first  look  at  the  condi¬ 
tions  required  to  satisfy  Eqs  (50)  and  (52)  ,  namely  the  elimina¬ 
tion  of  control  spillover.  At  the  core  of  this  method  will  be 


(50) 

(51) 


an  attempt  to  find  some  transformation  matrix,  T,  such  that 
subsequent  control  vector,  U(t),  required  for  Eq  (6)  will  be: 

U(t)  =  Tz(t)  (54) 

where  z(t)  is  now  the  new  control  input  and  with  constraint  that: 


B  T  =  0 
s 


(55) 


while: 


B  T  *  0 
c 


(56) 


One  method  with  which  to  obtain  this  transformation  matrix 
employs  a  technique  known  as  Singular  Value  Decomposition  (Ref  9) . 
Using  SVD  allows  reformulation  of  the  s  x  m  Bo  matrix  as: 

B  =  WZVT  (57) 

s 


where 


and 


W  is  an  s  x  s  orthogonal  matrix  of  left  singular 
vectors 

V  is  an  m  x  m  orthogonal  matrix  of  right  singular 
vectors 


S  '  0 

z  = .  (58) 

0  *  0 

_  _s  x  m 


Such  that  S  is  a  q  x  q  diagonal  matrix  of  the  singular  values  of 

B  ,  or:  (continued) 

s 
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Singular  values  are  always  greater  than  or  equal  to  zero, 
and  the  total  number  of  non-zero  singular  values  is  equal  to  the 
rank  of  the  decomposed  matrix.  As  long  as  Bg  is  full  rank  with 
dimensions  of  s  x  m,  then  q  is  the  minimum  value  of  the  pair 
(s,m) .  By  arbitrarily  letting  r  be  the  difference  between  q  and 
m,  or: 


q  +  r  =  s 


(59) 


The  W  can  be  partitioned  such  that: 


W  = 


Vf  :  V* 
q  .  r 


(60) 


S  X  s 


having: 


W  as  an  s  x  q  matrix 

q 

W  as  an  s  x  r  matrix 
r 


In  a  similar  fashion,  we  can  choose  p  as  the  difference  between 
q  and  n  such  that: 


q  +  p  =  m 


(61) 


We  can  now  partition  the  right  singular  vector  matrix,  V,  as: 

\ 


V  = 


v  :  v 
q  •  p 


(62) 


m  x  m 


V  as  an  m  x  q  matrix 

q 

V  as  an  m  x  p  matrix 

P 
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having : 


By  defining  as  our  transformation  matrix,  T,  we  find 
some  highly  desirable  results  with  respect  to  Eqs  (55)  and  (56) , 
namely: 


B  T  =  B  V  =  W  SV  V 
s  s  p  q  q  P 


(63) 


However,  since  V  is  an  orthogonal  matrix: 


V  V  =0 
q  P 


(64) 


Hence: 


BgT  =  0  (65) 

Coordinating  this  expression  for  the  transformation  matrix  with 

the  model  so  far  established,  it  should  first  be  noted  that  the 

dimensions  of  Bg  are  directly  the  result  of  both  the  number  of 

modes  to  be  supressed  (s)  and  the  number  of  actuators  employed 

(m) .  If  the  rank  of  the  matrix  Bg  is  equal  to  the  number  of 

actuators  available,  then  q  =  m.  By  Eq  (61)  it  is  seen  that  this 

forces  p  to  be  zero.  It  follows  that  V  =  T  =  0,  and  we  are 

P 

restricted  to  the  trivial  solution.  Recalling  the  previous  com¬ 
mittment  for  the  transformation  method,  this  would  fail  by 
allowing  BcT  =0.  It  is  clear  that  for  the  transformation  method 
to  be  carried  to  an  exact  solution,  special  conditions  in¬ 
cluding  q  <  m  must  be  met.  Restated,  the  rank  of  B  must  be 

s 

less  than  the  number  of  actuators.  It  should,  however;  be 
noted  that  if  you  are  restricted  to  a  fixed  number  of  pre¬ 
oriented  actuators,  performance  is  enhanced  by  using  the  sin¬ 
gular  vector  associated  with  the  least  singular  values  (even 
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though  B  T  #  0) .  In  any  case,  where  q  =  m,  an  order  reduction 
s 

scheme  is  required  to  get  an  exact  non-trivial  solution  to 

B  T  =  0.  As  a  verifying  set  of  examples,  let  q  be  the  rank  of 
s 

B  ,  and  let  that  matrix  have  dimensions  s  x  m  with  s  =  4  and 

5 

ra  =  3.  If  B  is  full  rank,  q  =  3.  Therefore  by  Eq  (61),  p  =  0. 
s 

However,  if  we  can  reduce  the  rank  of  B  to  q  =  2,  then  p  =  1 

s 

and  Vp  is  non-zero.  As  will  be  demonstrated,  this  rank  defi¬ 
ciency  is  obtained  either  through  judicious  orientation  of  the 
actuators  (driving  a  non-zero  singular  value  to  zero)  or  through 
addition  of  actuators  and  increasing  m.  The  minimum  number  of 
actuators  that  can  be  used  where  the  former  method  is  employed 
is  two,  since  a  matrix  of  rank  1  cannot  be  made  rank  deficient. 

Regardless  of  how  an  appropriate  non-zero  T  is  formed,  we 
will  how  have  the  resultant  solution  vector  in  Eq  (54)  ,  where: 


z  =  -GtXc 


Rfc  =  -T  RT 


such  that  Rt  is  a  p  x  p  positive  definite  matrix,  and: 


B  =  B  T 
t  c 


with  A  B.  completely  controllable. 

C  t  \ 

To  generate  the  new  control  vector,  z,  the  same  approach 
as  followed  in  Section  III  is  employed.  The  control  gain  matrix 
is  now  defined  by: 


-1  T 

Gt  "  Rt  Bt  Pc 
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and 


0 


(70) 


PA  +  A  TP  +  P  BJ.R71BJ.TP  +  Q  = 
cc  cc  ctttc  c 


where: 

is  a  p  x  m  reduced  degree  of  freedom  critical  state 
feedback  gain  matrix 

Pc  is  a  m  x  m  positive  definite  solution  to  the  reduced 
order  Ricatti  Equation 

The  gain  matrix  is  therefore  finally  transformed  by: 


Gc  =  TGfc  (71) 

which  will  produce  a  new  m  dimensional  control  with  zero  control 
spillover. 

A  parallel  technique  is  employed  to  eliminate  the  KCg 
observation  spillover  term.  Here,  the  number  of  sensors  must 
exceed  the  number  of  suppressed  modes,  or  Cg  must  be  made  rank 
deficient  through  sensor  re-orientation.  The  specific  methodo¬ 
logy  for  reducing  the  order  of  the  optimal  regulator  will  be 
described  as  part  of  the  computer  model,  and  in  the  investigation 
which  follows. 
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Computer  Model 

Appendix  A  represents  a  computer  listing  for  one  run  of  the 
main  program.  This  particular  run  applies  the  transformation 
method  to  an  eight  mode  nominal  model  such  that  the  control  spill¬ 
over  (BgG)  is  driven  to  zero.  Although  the  program  is  seen  to  be 
quite  lengthy,  the  comment  cards  which  have  been  included  for 
clarity  suggest  the  overall  straightforwardness  of  the  approach. 

As  a  first  step,  the  parameter  matrices  (A,  B,  C)  are  built. 
The  A  matrix  portion  of  the  program  reads  in  the  natural  fre¬ 
quencies  from  the  NASTRAN  data,  and  uses  these  frequencies  and 
a  prescribed  damping  ratio  (0.005)  to  fill  this  parameter  matrix 
appropriately.  The  B  and  C  matrices  are  formed  as  a  matrix  pro¬ 
duct  of  mode  shapes  and  actuator  or  sensor  locations.  That  is, 

B  is  formulated  as: 

B  =  4>TD  (72) 

where 

<)>  is  the  matrix  whose  columns  are  the  eigenvectors 
for  each  mode  supplied  by  NASTRAN  (mode  shapes) 

D  is  a  direction  cosine  matrix  for  the  locations 
and  orientations  of  the  prescribed  actuators. 

Since  the  sensors  and  actuators  are  collocated,  it  then  becomes 
clear  that: 

T  T 

C  =  B  =  D  $  (73) 

Next,  by  supplying  as  an  input  value  the  number  of  modes  to  be 
controlled,  the  program  takes  the  A,  B,  and  C  matrices  and 
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generates  their  controlled  and  suppressed  counterparts  (i.e.  Ac, 

A  , B  ...).  With  these  matrices  formed,  along  with  their  trans- 
poses,  the  steady  state  feedback  gain  matrices  (K  and  G)  are 
established.  Reviewing  this  process  as  described  in  Section  III, 
it  is  seen  that  one  step  involves  solution  of  the  steady  state 
matrix  Ricatti  Equation.  This  solution  is  obtained  via  highly 
specialized  computer  subroutines  created  by  Kleinman  (Ref  9) . 

With  these  gain  matrices,  the  total  system  equation  seen  as  Eq 
(44)  is  formed.  An  eigenvalue  analysis  using  subroutine  EIGRF 
from  the  International  Mathmatical  and  Statistical  Library 
(IMSL)  is  completed  against  the  controller  (A  +  BcG)  the  ob¬ 
server  (A  -  KCc) ,  and  the  entire  system.  This  allows  for  a 
stability  analysis  based  on  these  eigenvalues. 

Next,  a  time  history  response  (20  seconds)  is  performed  on 
the  line  of  sight  in  both  the  x  and  the  y  directions  at  grid  point 
1.  This  is  accomplished  in  two  steps.  First,  the  CC6600  sub¬ 
program  library  of  the  Air  Force  Institute  of  Technology  is 
implemented  such  that  program  ODE  (Ref  10)  can  be  used  to  inte¬ 
grate  the  state  equation: 

X  =  AX  +  Bu  (74) 

to  establish  x(t)  for  t  =  0.0,  0.1,  0.2, ...20.0.  Then,  using 
the  mode  shapes  and,  primarily  their  x  and  y  components  at  grid 
point  1,  the  line  of  sight  magnitudes  are  formulated  such  that: 

n 

X  (t)  =  Z  <f>  -x.  (t)  (75) 

1  *  _  ill 
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and  n 

X_(t)  =  E  <J>  .x.  (t)  (76) 

2  i=i  21  1 

with 

X^(t)  being  the  line  of  sight  in  the  x  direction 
X2(t)  being  the  line  of  sight  in  the  y  direction 

This  set  of  results  provides  a  baseline  for  comparison  of 

future  analyses.  Once  these  plots  are  completed,  a  singular 

value  decomposition  is  performed  on  Bg  as  the  first  step  in  the 

transformation  method.  The  actuator  corresponding  to  grid 

point  7  is  rotated  incrementally  until  the  least  singular  value 

of  B  becomes  nearly  zero.  In  effect,  this  reduces  the  rank 
s 

of  B  .  With  this  new  orientation,  a  new  control  gain  matrix 
s 

(G)  is  formed  using  the  methods  described  in  Section  III. 

Also,  new  B  and  C  matrices  are  created  to  account  for  the  re¬ 
oriented  sensor /actuator  pair. 

With  these  new  values,  the  program  returns  to  the  eigen¬ 
values  obtained  previously.  The  fact  that  the  system  eigenvalues 
are  those  of  the  diagonal  members  is  born  out.  New  plots  are 
then  generated  so  as  to  compare  the  time  history  responses  with 
and  without  control  spillover.  This  same  approach  is  followed 

in  driving  the  observation  spillover  (KC  )  to  zero.  This  set 

s 

of  runs  demonstrates  the  improvement  available  without  adding 
hardware . 

Finally,  two  sensors  are  added  at  grid  point  one  to  examine 
the  effectiveness  of  adding  some  fairly  simple  hardware  (as 
opposed  to  adding  actuators) .  These  two  sensors  are  given  an 
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orientation  typical  to  those  sensors  already  prescribed.  This 
run  is  repeated  against  a  twelve  mode  model  to  verify  the  legit¬ 
imacy  of  the  first  truncation  of  higher  modes.  It  should  be 
noted  that  the  selection  of  two  additional  sensors  at  grid  point 
one  was  arbitrary.  Any  number  of  additional  sensors  could  be 
added  at  any  location  for  this  final  study. 


Outline 


Investigation 


A  systematic  approach  toward  assessing  the  effectiveness 

of  the  transformation  method  was  initiated.  As  a  first  case, 

the  system  eigenvalues  and  line  of  sight  time  history  responses 

were  examined  for  models  with  and  without  control  spillover 

(B  G) .  The  transformation  technique  was  only  applied  to  the 

control  gain  (G) .  The  sensor/actuator  pairs  remained  collocated 

while  one  of  these  pairs  was  rotated  to  produce  an  additional 

zero  singular  value  to  Bg.  An  angular  orientation  was  obtained 

which  produced  a  rank  reduction  in  the  suppressed  control  matrix 

The  weighting  function  of  the  controlled  states,  Xc,  was  set  at 

the  identity  matrix.  Upon  successful  completion  of  this  first 

case,  the  process  was  repeated  with  increasingly  higher  control 

weighting.  Then,  this  set  of  runs  was  compared  to  the  case 

of  eliminating  observation  spillover  (KC  ) ,  rather  than  control 

s 

spillover.  The  purpose  of  this  alteration  to  the  main  program 
is  twofold.  First,  it  would  demonstrate  that  the  total  system 
matrix  (Eq  (44))  is  block  diagonalized  successfully  by  forcing 
either  of  the  spillover  terms  to  zero.  Second,  it  facilitates 
the  final  area  of  investigation;  namely  the  potential  benefits 
of  sensor  additions.  The  addition  of  sensors  (rather  than 
actuators)  within  the  prescribed  model  was  chosen  out  of  practi¬ 
cality.  From  a  "hardware"  viewpoint,  the  addition  of  sensors  is 
seen  to  be  considerably  more  realistic  than  the  addition  of 
point  force  actuators. 


For  all  cases  examined,  the  overall  attempt  is  to  reduce 
the  line  of  sight  error  in  the  x  and  y  directions  at  grid  point 
one  to  less  than  0.0004  radians  and  less  than  0.00025  radians 
respectively  in  20  seconds. 

Elimination  of  Control  Spillover 

In  an  attempt  to  further  clarify  the  direction  of  this 
analysis  recall  from  Sections  II  and  III  that  the  control  gain 
matrices  are  determined  using  steady  state  optimal  regulator 
theory,  which  involves  minimization  of  related  quadratic  perfor¬ 
mance  indices.  These  performance  indices  for  the  model  with  and 
without  spillover  are: 

J  =  *5/°°  (X?FX  +  uTRu)  dt 

0  c  C 

J  =  %f°°  (X  TFX  +  UR.  U)  dt 
0  C  Sr  t 

respectively. 

An  inspection  of  these  two  indices  demonstrates  the  role  of 
the  control  weighting  matrix,  F,  as  an  amplifier  of  the  resul¬ 
tant  gains  applied  to  the  controlled  states.  All  cases  run 
attempt  to  control  the  first  two  modes  and  suppress  the  remain¬ 
ing  six.  An  attempt  to  modify,  and  ultimately  improve  perform¬ 
ance  is  tied  to  increasing  the  magnitude  of  this  weighting 
matrix.  It  is  known  from  the  previously  developed  theory  that 
increasing  the  magnitude  of  the  control  gain  (here  G)  has  the 
coincident  negative  effect  of  increasing  control  spillover.  It 
is  with  this  awareness  that  the  first  study  is  accomplished. 

This  study  involves  generating  the  system  matrix  and  examining 
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the  eigenvalues  and  associated  line  of  sight  time  history  re¬ 
sponse.  Once  these  data  are  generated,  the  transformation 
technique  of  Section  III  is  applied  to  force  BgG  to  zero. 

Table  4  is  a  presentation  of  data  pertinent  to  the  first 
case  in  which  F  is  set  at  the  identity  matrix.  Both  sets  of 
system  eigenvalues  exhibit  stability.  Additionally,  eigen¬ 
values  of  the  entire  system  are  the  same  as  the  eigenvalues 
of  the  matrices  on  the  diagonal,  verifying  that  control  spill¬ 
over  has  been  eliminated.  It  should  also  be  noted  that  the 
transformation  method  has  generated  a  controller  (A  +  BcG  -  KC) 
that  is  unstable,  but  which,  none-the-less ,  produces  a  stable 
system.  Figures  6  through  9  represent  the  time  history  responses 
for  the  x  line  of  sight  and  y  line  of  sight  errors.  Although  a 
precise  bandwidth  on  the  error  is  difficult  to  establish,  it  is 
obvious  that  the  prescribed  limits  specified  in  the  outline 
portion  of  this  section  have  not  been  satisfied. 

The  next  step  then  involves  multiplying  the  F  matrix  by 
scalar  powers  of  ten  (i.e.  1,  10,  100,...).  Until  F  reaches 
1000[I],  there  is  no  significant  improvement  in  the  line  of 
sight  error  for  either  the  case  with  or  without  control  spill¬ 
over.  However,  at  the  F  1000[I],  significant  changes  in  the  sys¬ 
tem  response  become  evident.  Table  5  is  a  presentation  of  the 
associated  eigenvalues  for  this  case.  Clearly,  the  spillover 
terms  have  now  forced  the  system  (without  transformation  of  the 
control  gain)  unstable.  The  eigenvalues  after  gain  transforma¬ 
tion,  however,  still  exhibit  stability.  This  demonstrates  the 


certain  advantages  of  using  this  method.  Figures  10  and  11 
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depict  line  of  sight  errors  without  transformation  and  demon¬ 
strates  the  unstable  response  indicated  by  the  associated  eigen¬ 
values.  The  time  histories  of  Figures  12  and  13  present  an  x 

line  of  sight  error  for  the  system  with  B  G  =  0  within  an 

s 

approximate  bandwidth  of  ±.0013,  and  a  y  line  of  sight  error 
of  ±.0008.  As  the  control  gain  weighting  function  is  increasing 
there  is  no  significant  improvement  of  response.  The  trend  of 
these  data  suggests  that  the  criteria  for  pointing  accuracy 
cannot  be  met  with  the  prescribed  number  of  sensors  and  actua¬ 
tors  (6  each) .  It  is  clear,  for  the  reasons  highlighted,  that 
sensors  will  have  to  be  added. 


Table  IV 


Elimination  of  Control  Spillover;  F  =  l.oflj 


System  Eigenvalues 


Before 

Transformation  (B_G#0) 

s  ■ 


After 

Transformation  (B_G=0) 

s 


02822 

+ 

5 . 70935i 

-s 

-.02855 

+ 

5 . 71073i 

02838 

± 

5 . 67583i 

-s 

-.02838 

± 

5 . 67583i 

02553 

+ 

5 . 14848i 

-s 

-.02575 

+ 

5 . 14935i 

01918 

+ 

3 . 84804i 

-s 

-.01924 

± 

3 . 84834i 

01778 

+ 

3 . 55770i 

-s 

-.01778 

± 

3 . 55770i 

01467 

+ 

2 . 96372i 

-s 

-.01482 

± 

2. 96458i 

08663 

+ 

1 . 47902i 

-c 

-.07712 

± 

1 . 46602i 

06679 

+ 

1 . 18915i 

-c 

-.00751 

± 

1 . 17064i 

06279 

+ 

1 . 45703i 

-0 

-.08627 

± 

1 . 46583i 

03768 

+ 

1 . 16069i 

-0 

-.04420 

± 

1 . 17052i 

Eigenvalues  of 

V-o 

+  B  G 
c 

.07457 

+ 

1 . 46607i 

-c- 

-.07712 

± 

1 . 46602i 

.05199 

+ 

1 . 17046i 

-c- 

-.00751 

± 

1 . 17064i 

Eiqenvalues  of  A  -  KC 
— 3 - c - c 


-.07457  +  1 . 46607i 
-.05199  i  1 . 17046i 


-O-  -.08627  ±  1 . 46583i 

-O-  -.04420  ±  1 . 17052i 


Eiqenvalues  of  A_  +  B_G  -  KC_ 
c  c  c 


-.00733  ±  1 . 46222i 
-.00585  +  1 . 16818i 


.00194  ±  1 . 46119i 
.03072  ±  1 . 17082i 


C  =  Controlled  Mode  Eigenvalues 
S  =  Suppressed  Mode  Eigenvalues 
0  =  Observer  Mode  Eigenvalues 


Figure  6.  LOSX  VS.  TIME,  B  G  /  (J ,  F  =  1  0 

40  S 


Figure  7.  LOSY  VS.  TIME,  B  G  /  0,  F  =  1 
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CD 


TIME 

Figure  8.  LOSX  VS.  TIME,  BG=0,  F=1.0 
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Table  V 


Elimination  of  Control  Spillover;  F  =  1000.0  [ij 


Before 

Transformation  (B_G#0) 

s 

.09178  ±  5 . 73850i 
.02838  ±  5 . 67583i 
.03098  +  5 . 17892i 
.00835  ±  3 . 85687i 
.01764  ±  3 . 55783i 
.01328  ±  2 . 98115i 
.300048  +  Oi 
1.15711  +  Oi 
1.66482  ±  . 80887i 
.07067  ±  1 . 46398i 
.04480  ±  1 . 16631i 


System  Eigenvalues 


-S- 

-S- 

-s- 

-s- 

-s- 

-s- 

-c- 

-c- 

-c- 

-o- 

-o- 


Af  ter 

Transformation  (B  G=0) 
- — s - - 

-.02855  ±  5 . 71073i 

-.02838  +  5 . 67583i 

-.02575  +  5 . 14935i 

-.01924  ±  3 . 84834i 

-.01778  +  3 . 55770i 

-.01482  +  2 . 96458i 

-3.11793  +  Oi 

-1.46147  +  Oi 

-.02000  +  1 . 17188i 

-.08627  ±  1 . 46583i 

-.04420  ±  1 . 17064i 

-3.11793  +  Oi 
-1.46147  +  Oi 
-.02000  +  1 . 17188i 

-.08627  ±  1 . 46583i 
-.04420  ±  1 . 17064i 

-3.23401  +  Oi 
-1.18781  +  Oi 
-.01849  ±  1 . 17817i 


2.89759  +  Oi 
1.53062  +  Oi 
1.5091-  ±  . 74849i 


Eigenvalues  of  A  +  B_G 
c  c 

-c- 

-C- 

-c- 


-.07458  ±  1 . 46608i 
-.05199  ±  1 . 17046i 


Eigenvalues  of  A.  -  KCg 

-O- 

-O- 


Eiqenvalues  of  A  +B  G  -KC 
— a - c - c - c 


-3.02017  +  Oi 
-1.46304  ±  . 64937i 
-1.27355  +  Oi 


C  =  Controlled  Mode  Eigenvalues 
S  =  Suppressed  Mode  Eigenvalues 
O  =  Observer  Mode  Eigenvalues 


44 


LOSY  *10'z 

.-0.32  -0.24  -0.1B  -0.08  0.00  0.08  0.16  0.24  0.32 


0.00 


T 


- 1 - 1 - 1 - 1 - 

8.00  12.00  18.00  20.00 

TIME 

Figure  ll.  LOSY  VS.  TIME,  BgG  *  0,  F  *  1000.  U 
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±’.0015 


4.00 


24.00 
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Figure  13.  LOSY  VS.  TIME,  BgG  =  0,  F  ® 
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TIME 


±.0009 


24.00 


1000.0 


Elimination  of  Observation  Spillover 

The  approach  taken  during  this  portion  of  the  analysis  is 
directed  by  an  awareness,  a  priori,  that  sensors  will  be  added. 

As  a  preliminary  step,  the  procedure  for  applying  the  transforma¬ 
tion  method  to  the  control  gains  is  first  reapplied  to  the  obser¬ 
vation  gain  (K) .  The  same  sensor  and  actuator  pair  at  grid  point 
seven  is  again  rotated  until  the  smallest  singular  value  of  Cg  is 
driven  to  zero.  Table  VI  presents  the  results  of  the  eigenvalue 
analysis  which  followed.  Q  replaces  F  as  the  observation  weight¬ 
ing  matrix  acting  on  the  controlled  states.  Once  again,  the 
system  matrix  is  seen  to  be  stabilized  and  diagonalized  via  the 
transformation  method.  A  more  pertinent  case,  in  light  of  a  forth¬ 
coming  examination  of  sensor  additions,  is  an  application  of  the 
transformation  method  to  the  system  with  the  sensors  in  their 
original  orientation.  Table  VII,  below,  lists  the  singular  val¬ 
ues  of  C  for  the  fixed  six  sensors.  An  examination  of  their 
s 

relative  magnitudes  indicates  that,  although  the  last  singular 
value  is  non-zero,  some  potential  benefits  may  be  gained  by  apply¬ 
ing  the  transformation  method  with  this  singular  value  and  its 
associated  right  singular  vector. 

Table  VII 

Singular  Values  of  Cg;  Six  Non-reoriented  Sensors 


Number 

Singular  Value 

1. 

.70706 

2. 

.70423 

3. 

.70363 

4. 

.49803 

5. 

.42875 

6. 

.28536 
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Table  VI 


C  =  Controlled  Mode  Eigenvalues 
S  =  Suppressed  Mode  Eigenvalues 
0  =  Observer  Mode  Eigenvalues 
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Table  VIII 


Reduction  of  Observation  Spillover;  Q  =  1000. o£l]  (Fixed  Sensors) 


System  Eigenvalues 


Before 

Transformation  (KC_#0) 

s 

After 

Transformation  (KC_^ 
"  s 

.12469  ±  5 . 69223i 

”S“ 

-.06666 

+ 

5.71292i 

-.02838  ±  5 . 67583i 

-s- 

-.02838 

± 

5 . 67583i 

.06669  ±  5 . 15788i 

-s- 

.09387 

± 

5 . 21012i 

.00396  ±  3 . 85266i 

-s- 

-.05836 

± 

3 . 81382i 

-.01746  ±  3 . 55780i 

-s- 

-.01331 

± 

3 . 56096i 

.03719  ±  2 . 97542i 

-s- 

-.02278 

± 

2 . 96146i 

-.06713  ±  1 . 45779i 

-c- 

-.02024 

± 

1 . 47982i 

-.03931  ±  1 . 16093i 

-c- 

-.09759 

± 

1 . 10147i 

-2.33756  ±  . 32091i 

-0- 

-.27202 

± 

1 . 34862i 

-1.72715  ±  1 . 20427i 

-0- 

-.03520 

± 

1 . 20825i 

Eioenvalues  of  A_  + 

B  G 
— c— 

-.07457  ±  1 . 46607i 

-C- 

-.07457 

± 

1 . 46607 i 

-.05199  ±  1 . 17046i 

-c- 

-.05199 

± 

1 . 17046i 

Eiqenvalues  of  A_  - 

KC 

-2.89750  +  Oi 

-0- 

c 

-.20967 

± 

1 . 17020i 

-1.5308  +  Oi 

-0- 

-.04992 

± 

1 . 46369i 

-1.50921  ±  .  74854i 

-O- 

Eigenvalues  of  A.  +  B_G 

-  KC. 

1.45136  t  . 95906i 

c 

-.01187 

± 

1 . 48607i 

2.13224  ±  1 . 3163i 

.13171 

± 

1 . 18510i 

C  *  Controlled  Mode  Eigenvalues 
S  =  Suppressed  Mode  Eigenvalues 
0  =  Observer  Mode  Eigenvalues 
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LO$X  VS.  TIME 


Max 

Error:  +.0030 


0.00  4.00 


B.  00 


12.00  16.00  20.00  24.00 

TIME 


Figure  14.  LOSX  VS.  TIME,  KCg  f  0,  Q  =  1000.0 
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TIME 


±.0015 


24 .00 


f  0,  Q  =  1000.0 


S.  TIME 


The  results  of  this  analysis  bear  out  our  expectations. 

First,  some  observability  has  been  gained,  and  response  is  im¬ 
proved.  This  is  born  out  by  examining  Figs  14  through  17, 
which  are  the  line  of  sight  response  histories  with  and  without 
transformation.  However,  an  examination  of  Table  VIII  shows  that 
the  additional  observability  was  not  sufficient  to  generate  a 
completely  stable  system.  The  observation  spillover  has  been 
reduced,  but  not  eliminated. 

Sensor  Additions 

Until  now,  we  have  seen  that  the  transformation  method  can 
be  successfully  applied  to  the  tetrahedron.  Spillover  can  be 
minimized  or  completely  eliminated,  depending  on  whether  or  not 
sensor  and  actuator  reorientations  are  permitted.  In  reality, 
it  is  perhaps  more  likely  that  one  would  have  less  than  complete 
liberty  to  do  this.  Regardless,  the  specified  line  of  sight 
criteria  has  not  been  met.  Hence,  we  are  left  with  sensor  addi¬ 
tions  as  a  last  resort. 

Two  sensors  were  added  to  the  system  at  grid  point  one. 

This  number  and  location  are  essentially  arbitrary,  but  will 
serve  as  a  starting  point  for  more  exhaustive  subsequent  analy¬ 
ses.  Table  IX,  using  the  format  applied  throughout  this  report, 
presents  the  results  of  this  case.  It  is  clear  that  sensor  addi¬ 
tions  have  allowed  the  same  system  matrix  block  diagonalization 
as  previous  techniques.  However,  as  has  been  the  case  previously, 
Figs  18  through  21  demonstrate  that  the  criteria  for  line  of 
sight  response  has  not  been  met.  The  improvement  to  note  is 
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5 

I 

between  time  history  response  associated  with  Table  VIII  and 
that  of  Table  IX.  Clearly,  the  addition  of  sensors  has  enhanced 
the  overall  performance. 


Table  IX 


Elimination  of  Observation  Spillover;  Q  =  lOOO.O^l"] 

System  Eigenvalues 

Before  After 


8  Sensors 


Transformation  (KCg^0) 


Transformation  (KCs=0) 


.01473  ± 

5 . 71307i 

-s- 

-.02855 

+ 

5 . 71071i 

.02838  ± 

5 . 67583i 

-s- 

-.02838 

± 

5 . 67583i 

.02350  ± 

5 . 15484i 

-s- 

-.02575 

± 

5 . 14935i 

.05240  ± 

3 . 84096i 

-s- 

-.01924 

± 

3 . 84834i 

.01538  ± 

3 . 55804i 

-s- 

-.01779 

± 

3 . 55770i 

.02363  ± 

2.9613U 

-s- 

-.01482 

± 

2 . 96458i 

.07371  ± 

1.46483i 

-o 

-.07457 

± 

1.466071 

.05187  ± 

1 . 17031i 

-o 

-.05199 

± 

1 . 17046i 

15.6899 

+  Oi 

-O- 

-.03632 

± 

1.193061 

5.60922 

4-  Oi 

-O- 

-.00733 

± 

1 . 46676i 

1.20084  ■ 

+  Oi 

-0- 

1.02196 

f  Oi 

-0- 

Eigenvalues  of  A  +B  G 

C-  w 

.07457  ± 

1 . 46607i 

-o 

-.07457 

± 

1 . 46607i 

.05199  + 

1 . 17046i 

-o 

-.05199 

± 

1 . 17046i 

Eigenvalues  of  A_  -  KC_ 

15.69921 

+  Oi 

-  -  1  **■  -  .  c  c 

-0- 

-.03632 

± 

1. 19306i 

5.62028 

+  Oi 

-0- 

-.00733 

± 

1 . 46676i 

1.15106 

+  Oi 

-0- 

1.01144 

+  Oi 

-0- 

Eigenvalues  of  Ac  +B_G  - 

KC 
— c 

15.33771 

+  Oi 

-.07457 

± 

1 . 46231i 

4.12626 

+  Oi 

-.02153 

± 

1 . 19748i 

2.48129  +  Oi 
1.25725  +  Oi 
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Figure  18.  LOSX  VS.  TIME,  KCs  *0,8  Sensors,  Q  =  1000.0 
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Figure  20.  LOSX  VS.  TIME,  KC  =  0,  8  Sensors, 
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Figure  21.  LOSY  VS.  TIME,  KCg  =  0, 
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Conclusions 


Two  key  conclusions  can  be  drawn  from  the  preceding  analy¬ 
ses.  First,  given  a  fixed  number  of  sensors  and  actuators  with 
fixed  orientation,  the  destabilizing  effect  of  observation  spill¬ 
over  and  control  spillover  can  be  "minimized".  When  a  reorienta¬ 
tion  of  those  sensors  and  actuators  is  permitted,  these  spillover 
effects  can  be  completely  eliminated.  Elimination  of  either  con¬ 
trol  spillover  or  observation  spillover  guarantees  system  stabil¬ 
ity,  regardless  of  whether  or  not  response  criteria  are  satis¬ 
fied.  Second,  if  sensor  reorientation  is  not  allowed,  complete 
elimination  of  observation  spillover  can  still  be  accomplished 
through  sensor  additions. 

The  transformation  method  was  found  to  be  very  effective  in 
eliminating  control  spillover  and  uncoupling  system  eigenvalues 
when  the  number  of  actuators  in  the  system  is  greater  than  the 
number  of  modes  to  be  suppressed.  When  the  number  of  modes  to 
be  suppressed  is  equal  to  the  number  of  actuators,  complete 
elimination  of  control  spillover  can  be  accomplished  through 
an  actuator  reorientation  which  reduces  the  rank  of  the  control 
matrix,  B.  A  parallel  case  can  be  made  for  the  elimination  of 
observation  spillover  where  the  number  of  sensors  is  greater 
than  or  equal  to  the  number  of  suppressed  modes.  When  reorien¬ 
tation  is  not  permitted,  the  degree  of  response  improvement  is 
strictly  a  function  of  the  relative  magnitudes  of  the  singular 
values  of  the  decomposed  matrices.  For  the  specific  cases  exam¬ 
ined,  the  truncation  of  higher  frequency  modes  was  seen  to  be 
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Recommendations 


The  major  theme  of  this  analysis  suggests  that,  due  to  the 
complexity  of  larger  and  larger  space  systems,  controllers  will 
have  to  be  developed  to  operate  on  only  those  modes  critical  to 
system  response.  This  requirement  is  imposed  due  to  limited 
computer  and  hardware  capabilities.  Since  line  of  sight 
was  established  as  the  performance  criteria  in  this  study,  the 
modes  were  arranged  in  order  of  decreasing  displacement  at  the 
selected  grid  point.  The  decision  to  control  two  modes  and  to 
suppress  six  was  arbitrary.  Since  the  selection  of  "critical 
modes"  is  the  starting  point  in  developing  an  eventual  control¬ 
ler,  the  importance  of  this  step  cannot  be  overemphasized.  No 
automated  technique  for  this  process  is  currently  available.  An 
exhaustive  re-application  of  the  computer  technique  found  in 
Appendix  A  may  result  in  satisfaction  of  the  prescribed  time  re¬ 
sponse  criteria.  More  importantly,  valuable  insight  into  this 
task  of  critical  mode  selection  might  be  obtained  as  fallout 
from  this  study.  In  a  parallel  sense,  the  selection  of  two 
sensors  to  be  added  for  the  final  case  examined  was  also  ar¬ 
bitrary.  Once  again,  a  follow  up  with  varying  numbers  and  loca¬ 
tions  of  additional  sensors  would  be  necessary  to  develop  the 
optimal  controller  for  this  model.  Finally,  sensitivity  to 
modelling  inaccuracies  would  be  a  natural  topic  for  further 
analysis.  Parameter  variations  would  have  to  be  incorporated 
into  the  NASTRAN  analysis  provided  in  order  to  simulate  mode 
shape  and  frequency  errors  for  this  sensitivity  study. 
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Appendix  A 


Eigenvector  Results  of  NASTRAN  Analysis 


Eigenvalue  1  = 

Eigenvector  1  = 


Eigenvalue  3  = 

Eigenvector  3  = 


Eigenvalue  5  = 

Eigenvector  5  = 


Real  Eigenvectors 


1 . 37043E+00 

Eigenvalue  2  - 

2 . 15145E+00 

-2 . 47073E-01 
4 . 27857E-02 
1 . 45180E-06 
-1 . 96263E-02 
3 . 39753E-02 
-7 . 21326E-02 
-3 . 69602E-02 
4 . 39747E-02 
-1.96224E-02 
5 . 29624E-02 
_4 . 39672E-02. 

Eigenvector  2  = 

"  3 . 99896E-01 

2 • 30929E-01 
-1.48908E-01 

8 . 32862E-02 

4 . 80849E-02 

6 . 81283E-02 

6 . 99996E-02 

2 . 25294E-02 
-4 . 72104E-02 
5. 45051E-02 

4 . 93610E-02 
-4 . 72153E-02_ 

8. 78894E+00 

Eigenvalue  4  = 

1 . 26576E+00 

"*6 . 36794E-02 
3 . 67778E-02 
4.00015E-01 

1 . 98377E-01 

1 . 14530E-01 
2 . 00976E-01 

1 . 54760E-01 

6 . 80356E-02 

9 . 78233E-02 

1 . 36292E-01 

1 . 00014E-01 

9 . 78391E-02 

Eigenvector  4  = 

"  2.74559E-02 
-4 . 75782E-02 
-2 . 24915E-05 
-1 . 71840E-01 

2 . 97744E-01 
-6 . 81682E-05 
-2 . 51248E-01 
3.43581E-01 
-8 . 19014E-02 
-1.71848E-01 

3 . 89435E-01 
8.1921 8E-02 

1 . 4  8101E+01 

Eigenvalue  6  = 

2 . 65165E+01 

-8 . 78330E-02 
-5 . 07014E-02 
-1 . 29877E-01 

3 . 09503E-01 

1 . 78636E-01 
-3 . 51419E-01 

2 . 86593E-01 

1 . 22432E-01 

1 . 13906E-02 

2 . 49398E-01 

1 . 86846E-01 

1 . 14008E-02 

Eigenvector  6  = 

f 1.35323E-05_ 
.1 -2I8I6E-11 

3 . 40156E-11 
-2 . 04139E-01 

3 . 53548E-01 
-6 . 05706E-06 
-2.04139E-01 
-3.53548E-01 

1 . 08602E-04 

4 . 08202E-01 

6 . 80214E-10 

5 . 06531E-10 
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Eigenvalue  7 
Eigenvector  7 


Eigenvalue  9 
Eigenvector  9 


Eigenvalue  11  = 
Eigenvector  11= 


Eigenvector  =  lx 


3 . 22159E+01 

Eigenvalue  8  = 

3 . 26133E+01 

-2 . 66140E-02~ 
4.60655E-02 
3 . 30215E-05 
3 . 37411E-02 
-5 . 84417E-02 
3.23144E-05 
2.73330E-02 
-5. 48104E-02 
-4 . 91269E-01 
3 . 38171E-02 
-5 • 10814E-02 

Eigenvector  8  = 

-2.99367E-02* 
-1 . 73Q93E-02 

8 . 78423E-02 

4 . 07052E-02 

2 . 35996E-02 
3.55373E-02 
2.74211E-02 

2 . 79794E-02 
-4 . 87453E-01 

3 .  79914E-02 

9 . 80954E-03 

.  4.90852E-01_ 

-4.87867E-01 

7 . 99170E+01 

Eigenvalue  10  = 

1.06164E+02 

"  9.90668E-02" 
5 . 72029E-02 
1 . 72892E-01 
1 . 07566E-01 
6 . 21328E-02 
-4 . 95312E-01 
-1 . 67880E-01 
-2. 19818E-01 
-1 . 11010E-02 
-2.74347E-01 
-3 . 55381E-02 

Eigenvector  10= 

-3.38986E-03 

5 . 84999E-03 
-1.60534E-05 
-2. 28617E-01 
3.95968E-01 
4.96376E-05 
3.78349E-01 

4 . 55436E-02 

1  -1 . 47053E-02 
-2. 28600E-01 
-3.04859E-01 

L-1.10861E-02_ 

1.47172E-02 

1 . 19320E+02 

Eigenvalue  12  = 

1.95068E+02 

6 . 36959E-02 
3.67781E-02 

9 . 58836E-02 
-2 . 40062E-01 
-1.38592E-01 
-2.60496E-01 
-8. 60592E-02 

3 . 94412E-01 

6 . 96952E-03 

2 . 98410E-01 
-2 . 71939E-01 

Eigenvector  12= 

’  3.20580E-02’ 
1.85105E-02 

6 . 43806E-02 
-4.02579E-01 
-2 . 32435E-01 
-1 . 30450E-01 
3.20382E-01 
-1 . 58741E-01 
-9.27787E-03 

2 . 27168E-02 
3.56828E-01 

6 . 97Q73E-03 

-9 . 28169E-03 

1  >y1  'zi  >  ■  •  •  . 
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Appendix  b 

Main  Program  Listing 
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PROGRAM  TEfRA 


74/74  OPT  =  1 


FTN  4»  8*5i8 


PROGRAM  TETRA  (1NPUT= /8G  ,  OtJT=>UT=/ 1 32 ,  T  APE  5 ,  7  A  PE6  ,  T  A  PE  7) 

DIMENSION  Yl2C),YP(2  0),WORK(t>2C)  ,  I  WORK  (51 
DIMENSION  CCN8(8,lb)  ,CS  P8  (  e,  16)  ,  CCNT  ( 16  ,  B >  ,  CT  (1  6  ,  1 ) 

DIMENSION  T1(1,6),TT1(6,1),V2(15,L),:TT(1,15> 

DIMENSION  VV2(16,16)  ,S2(ie,16),U.<2U6,15>  }uKl(l,16) 

DIMENSION  0K2(1, ib) ,GKT (16,16) 

DIMENSION  XXI <20 3) ,XX2(203> ,XX3(2D3)  ,  XX4(2t;‘) ,XXE(2C  3) ,  XX5 (203) 
DIMENSION  XX7  (203)  ,  X  X8(  23  3) ,  TM(  2:  3)  ,  X  1  <  2  J  3  i ,  X2  (  2C  3 ) 

DIMENSION  76 ( 2U, 20) 

DIMENSION  T  2  ( 2 , 8 )  , TT 2(8,2), 3T2(15»2)  ,  RT  2(2,8), RT3(2,  2) 

DIMENSION  V  3  ( 4 , 2 )  ,CTT2(2,16)  ,  VV3  (  ♦,  '♦)  ,  S3  (lb,  16)  >  ULX3  (16, IS) 
DIMENSION  RTII(2, 2>,OK3(2,lb)  ,D<4(2,16)  ,  3<T2(8,4> 

DIMENSION  B0(it»16)»SV(8)  ,M_ ( 16) 

DIMENSION  R 1  ( 6  ,  $ ) ,RK1  (8 , 8) ,  RRRl  ( 8 , 3 ) 

DIMENSION  PHI  (12,16)  ,PHIT(15,12)  ,3(12, 6), 31(6, 12) 

DIMENSION  PHIT5  ( 6 , 12 )  ,PHI5(12,6)  , B< { 3 , 6) ,  BnG  (8,  4)  ,  CR  (8 , 8) 
DIMENSION  KCR( 16,8) 

DIMENSION  CSP4 (12 ,8) 

DIMENSION  AT0T<16,16) ,BTOT(i8,6) ,:ror (fl, 15) 

DIMENSION  AC0N(16,16)  ,  AS  UP  (16, 13)  ,  ICON  (lb  ,L>  ,8SUF(  16  ,6) 
DIMENSION  CSP9 ( 16  ,6) 

DIMENSION  CY(16,16),ACN6(16,16) 

DIMENSION  WK (32) , CCT ( 16 , 16) , CN ( 15 , l 6) 

DIMENSION  WK1(32) , WK 2 (3 2) , W< 3 (32) 

DIMENSION  ACN1(16,16) ,UC(16,16) 

DIMENSION  CCON<8,lb)  ,CSUP(6,16) 

DIMENSION  CCN2 ( 8 , 16) ,CCN1  (15 , 8) , 83 N1 ( 16, b) 

DIMENSION  ACN2(lb,16)  ,CX  (16,  16)  ,1)LX(16,15) 

DIMENSION  3CN2(6, 16)  ,CSP1(8,16)  ,83* 1(16, 5) 

DIMENSION  CGM3 (16,16) ,0(16,16) 

DIMENSION  R (6, 6), 0(16, 16), RY (6, 15),  A30NT( 16, 16) , CC0NT(16,3) 
DIMENSION  RX(6,16) , P ( 16 , 16) ,0KTRN(3,l6) ,DK(lt,5)  ,RR(6,6) 

01 MENSION  BTT(l,16),UL(l6,15)  ,CC(l:>,lb)  , OS  PI LL ( 16, 16 ) 

DIMENSION  CSP1LL  (lb,  :’.6>  ,  rD(t  6,16) 

DIMENSION  C0M1(16,16/ , COM 2(16,15)  ,3(16,15) , GO (6, 16) 

DIMENSION  00(16,16), COC (16,1b) 

DIMENSION  SS(16,16), RRR (6,6) ,RU(15,8) ,RV(lb,6) 

DIMENSION  77(15,16), S0(l6, 15)  jBCOMT ( 5 , 15 > 

DIMENSION  PHIT1(16,12),8SP3(16,5) ,BSP4(15,e) ,PHTT2 (16, 12) 

01  MENSION  T(6,l)  ,TT(1,6)  ,  3T  (  16,1)  ,Ui(l,  5) 

DIMENSION  RT (1,1) , VI (16,1), VV (16,1 5)  ,  Si (16 ,1b) 

DIMENSION  ULX1 ( 16,16) ,QQ (15,16), RTI (1,1) ,G01 (1,15) ,G02 (1,16) 
COMPLEX  W  (2C )  ,70(20, 2o)  ,  WW  (3  )  ,700(8,8) 

COMPL  EX  Wl(b) ,W2(8),W3(8)  ,71 ( 8, 8 ) , 72 (  8 , 6 ) ,  73  (6,8) 

INTEG  ER  FFF,Fi,F2,F3,F4,  GGG,  G1,G2,G3,  G»* 

INTEGFR  N5,Nb,N7 

INTEGER  K7 

INTEGER  NB,KK 

INTEGER  F5,N7 

INTEGER  IFLAG, NEON , J J 

INTEGER  NY 

INTEGFR  INI T 

INTEGER  I Z,NN,I JOB 

INTEGER  L , MO, NO, I A, I 0,1 C 

INTEGER  I,J,N,M,F,G,FF,GG 

INTEGER  NDIM,N0IM1,KIN,  KOUT,  KPUNCM 

REAL  TI,TOUT , RELERR, ABSERR 


PROGRAM  TETRA 


74/7R  0PT  =  1 


PTN  4.8+516 


REAL  Cl, Q2, 03, 04, 05 

REAL  06 

REAL  TOL 

REAL  CAMF 

REAL  FO 

COMMON  Z9  (2L  ,2t  ) 

COMMON/MA IM/NDI M ,  NO  I  Ml ,  C0M1 

COMMON/INOUT/KOUT , KIN, K PUNCH 

COMM ON/MAIN 2/COM  2 

COMMO  N/MAIN3/C0M3 

EXTERNAL  XOOT 

READ*  ,N5 

NOIM= 16 

N0IM1 =17 

TOL  =  .  om 

KIN=5 

K0'JT=6 

KPUNC  H=7 

IER=0 

0 AMP=  .005 

NY  =  0 

N7  =  5 

KK  =  0 

M=8-N5 

F5=2*  N5 

S  =  2+M 

EE =N5  +1 

GG  =  M  +  1 

F1=F5  +1 

F2  =  2+  F5 

F3=F2  +1 

F4=F2+G 

IF{<K  .GT.OGO  TO  37 
READ* , Nb 
READ* ,N7 


CREAT  E  A  MATRIX 


00  20  1*1, it 
00  10  J*l,16 
ATOT ( I,J)=J.O 
CONTINUE 
CONTI NUE 
00  3u  1*1,8 
J=I+8 

ATOT ( I, J) *1.0 
CONTINUE 
00  35  1*3,16 
J*I-8 

REAO* , ATOT(I , J) 

CONTI NUE 
00  36  1*9,16 
REAO“  ,ATOT (1,1) 

ATOT (I, I) *-0AMP*2 .*A  TOT (1,1) 
CONTINUE 


72 


12-; 


PROGRAM  TETRA  74/74  0PT=1 


FT  N  4.6+518 


37  PRINT *,"  " 

PRINT*,"  " 

PRINT*,”  " 

PRINT*, "THIS  RUN  REPRESENTS  AN  ANA.YSIS  c3R  AN  EIGHT  NODE" 
PRINT*, "APPROXIMATION  TO  THE  SYSTEM, WITH  ** ,  N  5 , "  MOOES" 
PRINT  VCONTROLLEC  ANO  ",M,“  MOOES  SUPPRESSED" 

PRINT*,”  •* 

IF(N7.GT.O) GO  TO  169 
PRINT*,"  “ 


126 


CREATE  B  MATRIX 


130 


135 


140 


145 


161 


155 


161 


165 


no  so  j=i,b 

9T0T  (  I,J)=li.(/ 

50  CONTINUE 
7  1  CONTINUE 

00  72  1=1,16 
DO  71  J=i , 16 
C(I,J)=0.0 

IF (I • EQ.J)C(I,J)=1.J 
7i  CONTINUE 

7  2  CONTINUE 

DO  74  1=1,12 
00  73  J=1,E 
0(1, J)=Q.O 
7  3  CONTINUE 

74  CONTINUE 

DO  76  1=4,6 
00  75  J=1 , 2 

REAO* ,0(I,J) 

75  CONTINUE 

75  CONTINUE 

DO  78  1=7,9 
DO  77  J=3,4 
READ*  ,0(1 ,J) 

77  CONTINUE 

78  CONTINUE 

00  8C  1=16,12 
DO  79  J=6,6 
READ*  ,0(1, J> 

79  CONTINUE 

80  CONTINUE 

00  82  1=1,12 
DO  91  J=1 ,8 
PHI  (I  ,J)=0.(> 

51  CONTINUE 

82  CONTINUE 

00  94  J=9, 16 
00  83  1=1,12 
REAO*  ,PHI (I, J) 

83  CONTINUE 

84  CONTINUE 

00  86  1=1,16 
DO  85  J=1 , 1 2 

7  3 


17C 
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175 


18L 


PHIT( I,J)=PHI(J,I> 

83  CONTINUE 

35  CONTINUE 

37  L 1  =  1  f  • 

MO  =  i  2 
NO=5 
I A  =15 
I  0  =  12 
10  =  15 

CALL  VMULFF(PHIT,0,L1,M0,N0,IA,I3,3TDT,IC.  liIR) 


185 


19k' 


195 


200 


205 


2  HI 


215 


2  2 


0  CREATE  C  MATRIX 

C 

IF(NY  .GT.DGO  TO  58 
121  IF (NS  »GT. b) GO  TO  113 
no  111  1=1, b 
00  110  J=l,  12 
01 (I, J)=D( J,I> 

110  CONTINUE 
til  CONTINUE 
118  L 1  =  N  8 

MO  =  1 2 
N0  =  16 
I A  =  8 
19  =  12 
IC  =  3 

CALL  VMULFf-  ID1, PHI, LI, MO, NO,  IA, IB, CTOT ,1 C , IE R> 
00  li,  6  1=1,  NS 
DO  145  J=l,  8 
CTOT ( I,J) =CTOT (I , 8+J ) 

14  5  CONTINUE 
14F  CONTINUE 

DO  14  8  1  =  1,  NS 
00  147  J=9» 16 
CTOT  (  I,J)=E.C 
14^  CONTINUE 
148  CONTINUE 
GO  TO  150 

113  DO  115  1  =  1,2 
00  114  J=l, 12 
01 (I, J) =o ♦ l 

114  CONTINUE 
1 1 F  CONTINUE 

01 (1,  1)  =  ,  3535573906 
01(1, 2)=-. 6122637618 
01(1, 3)=,7C7lce7612 
01(2, 1)=-. 3535533906 
01(2, 2)  =  *  61 22637618 
01(2,  3)  =•  7 j7 lt67 812 
00  117  1  =  3,8 
00  116  J  =  1,  12 
01 (I, J)=0(J,I-2) 

116  CONTINUE 

117  CONTINUE 
L 1  *N6 


4 


. 8+513 


74 


PROGRAM  TE^RA  74/ 74  OPT  =  1 

MO  =  l? 

231  N0  =  16 

I  A  =  9 
19  =  12 
IC  =  9 

CALL  VMUIFF(D1,PHI,L1,M0,N0,IA,16,CTQT , 
2  3';  no  94  1=1,8 

DO  53  J=1 ,o 

CTOT ( I« J) =CTOT (I , J+d) 


53 

CONTI NUE 

3  4 

CONTINUE 

24 1: 

00  3b  1=1,6 

00  55  J=9,ll 

CTOT ( I,J)=O.C 

35 

CONTINUE 

3  6 

CONTI NUE 

245 

150 

j 

IF(N7.GT.C)G0  TO  134 

z 

•> 

CREATE  A  CONTROLLED  AND  A  SUPPRESS 

25fi 

M 

*> 

S 

169 

00  180  1  =  1,  N5 

00  17!)  J=1,N5 

AC0N(I,J)=AT0T(I,J) 

17  * 

CONTI NUE 

253 

180 

CONTINUE 

oo  ?nn  i=i ;N6 

00  19  0  J=FF,F5 

ACON  (  I,J)=ATOT(I,H*J) 

190 

CONTI NUE 

26f 

200 

CONTI NUE 

DO  229  I  *  F  F  »  F  6 

DO  219  J=  1 ,  N 5 

AC0MH,J)=A10T(I+K,J> 

21 P 

CONTINUE 

265 

22  o 

CONTINUE 

00  240  1=  FF , F5 

DO  230  J=FF,F5 

ACONt  I,J)=ATOT(I*M, J+M) 

230 

CONTINUE 

27C 

240 

CONTI NUE 

54 

DO  260  1=1, M 

00  25  0  J*1,N 

ASUP(I,J)=ATOT(I*Nt»,  J*N5> 

25  0 

CONTI NUE 

275 

26  0 

CONTINUE 

00  28  5  1  =  1, M 

00  27  0  J=GG, G 

ASUP ( I  ,  J) =ATOT (I +  N5* J*2*N5) 

27  o 

CONTINUE 

280 

28  0 

CONTINUE 

00  3C  9  1  =  GG , 6 

00  29  8  J*  1,  H 

ASUPCI,  J>=AT0T<2*N5*I,N5  +  J> 

29  P 

CONTINUE 

285 

330 

CONTINUE  75 
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DO  32  0  I  =  GGf G 
no  310  J=GG,G 

ASUP ( I,J)=ATOT  <2*Nb*I,2*N5+J) 
310  CONTINUE 
291'  32"  CONTINUE 


CREATE  0  CONTROLLED  AND  B  SUPPRESSED 


295 


3  JC 


3  3b 


310 


31b 


33  DO  349  1=1, N5 
00  330  J=1 , 6 
BCON(I,J)=BTOT(I,J) 

330  CONTINUE 
340  CONTINUE 

00  3b  0  I  =  FF,F5 
00  3*>  D  J=  1, 6 
0CON(I,J>=BTOT<I*M,J) 

35  "  CONTINUE 
350  CONTINUE 
103  00  38  0  1  =  1,  M 

DO  37  C  J=1,F. 
nSUP(I,J)=9T0T<N5*l, J) 

37  0  CONTINUE 
390  CONTINUE 

00  4G9  I  =  GG , G 
00  390  J=l,6 
9SUP(I, J)=BT0T(2*N5*1, J) 
390  CONTINUE 
490  CONTINUE 

IFtNY.GT.DGC  TO  561 


32C. 


CREATE  C  CONTROLLED  AND  C  SJPPRESSEO 


325 


330 


335 


340 


104  DO  V20  1=1, N6 
00  413  J=1,N5 
CCON ( I, J) =CTOT ( I , J) 

41 P  CONTINUE 
420  CONTINUE 

DO  44  0  1  =  1, N6 
DO  430  J=FF,F5 
CCON  < I, J) =0 • 3 
♦30  CONTINUE 

44  0  CONTINUE 

162  00  460  1  =  1, N6 

00  450  J*i,H 
CSUP < I, J)=CTOT(I, N6*J> 

45  0  CONTINUE 
460  CONTINUE 

00  480  1=1, N6 
00  47  C  J=GG,G 
CSUP( I,J) =CT  OT (I, 2*N5*J) 
47  0  CONTINUE 
480  CONTINUE 


76 


PROGRAM  TETRA 


74/74  OPT  =  1 


FTN  4. 6*518 


CREATE  WEIGHTING  MATRICES 


r* 


165 


49  9 
491 


49? 
49  3 

5  31 


49  6 
497 

63  2 


51  0 
520 
522 


5  3  0 

540 

541 


55  0 

560 

561 


57  0 
59  0 


30  0 
601 


IF(NT  »GT •  o>  GO  TO  541 
00  491  1=1*6 
00  49  0  J=l»6 
R(I,J)=U.O 

IP(I . EQ.  J) R<I* J) =1.0 

CONTI NUE 

CONTINUE 

IF  (N6  .GT.  6)  GO  TO  5ci 
90  =  1.  “ 

DO  49  3  1  =  1*6 

00  49  2  J=l*t 

R1  (I »  J)  =t  *  U 

IF (I . EO. J) R1 (I, J) =1. S 

CONTI NUE 

CONTI NUE 

GO  TO  502 

DO  497  1=1,8 

00  49  €  J  =  i,8 

R1 (I, J) *C • 0 

IFCI. EO.J) R1(I, J) =1.0 


CONTI NUE 
CONTINUE 
39  =  1.  — 

00  520  1=1, F5 
00  510  J=1,FS> 

0(1, J  )  =  0.C 

IFCI.  EQ.  J)  0  (I ,  J)  =09 

CONTINUE 

CONTINUE 

DO  54  0  1=1, F5 

00  530  J=1,F5 

ACONT  (I, J) =f CON (  J  ,1' 

CONTINUE 

CONTINUE 

DO  56  0 

00  550  J=1,N6 
CCONT  (I, J) =CCON ( J , 1) 
CONTI NUE 
CONTT  NUE 

IF (N7  .GT . 6) GO  TO  625 

00  58  fl  1=1,6 

00  570  J=  1,  F5 

BCONT  (I, J) =BC0N ( J ,1) 

CONTINUE 

CONTINUE 

PRINT*, ‘•GO** 

IF (NY  .GT . 10 GO  TO  562 


00  50 1  1*1,6 
00  SQO  J*l, 6 
RR (I , 
CONTINUE 
CONTINUE 
00  610  1*1,6 


i ; 

j  ! 

i: 

u 


77 


P*OG*AM  TEtRA 


7 4/74  OPTsi 


“TN  4, 8*318 


4a-! 


4  06 


41J- 


41? 


4  2L 


4  2f> 


43i» 


435 


R'?(I,  I)*-1./R(I,I) 

51  n  r  ONTINUE 

00  5 1 2  1=1, e 
00  511  J=  1 , 6 
RRR(I  ,J)=-Kf>  Cl, J I 
511  CONTINUE 
51?  CONTINUE 

I^tNe  .GT.6)G0  TO  619 
00  51  A  1=1, t 
DO  51 3  J=1 ,E 
fiRKJ  tJ)  =  6'<J 
517  CONTINUE 
61 A  CONTINUE 

00  515  1  =  1,  h 
RR1(J  ,I)=-1./RU,  I) 
61?  CONTINUE 

00  518  1=1, e 
DO  61 7  J=l,f 
RRR1( I,J)=-Fk1(I, J) 

517  CONTINUE 

518  CONTINUE 
GO  TO  625 

619  00  621  1=1,6 
DO  S2f)  J=l,6 
RRl(T,J)=0,r 

620  CONTINUE 

621  CONTINUE 

DO  622  1*1,8 
RR1<T,I>=-1./R1U,I) 

52  2  CONTINUE 

DO  62 A  1=1,8 
DO  62  3  J=l, 8 
RRR1 ( I, J) =-RRl (I, J) 
52  3  CONTINUE 

62*  CONTINUE 


SOL /E  RESPECTIVE  RIC.'jfTI  iOJATIOVS 


AA'J 


A  46 


At iO 


455 


525  ll=F5 
M0=N6 
N0=N6 
I A =16 
19=9 
IC  =  16 

CALL  VMULFF (CCONT ,RRR1, LI , HO , NO, 1 1, 13 ,RU, iC, 1ER) 
N0=F5 

CALL  VMULFF <RU, CCON, LI, HO, NO, I  A,  IB, S3, IC, TER) 
ir(N*  .GT.WGO  TO  637 

526  DO  52  8  1=1, F5 
00  527  J*1,F5 

ACN1 ( I, J) =ACON (I , J) 

ACN2 ( I, J) =ACONT (I , J) 

827  CONTINUE 
628  CONTINUE 
582  00  530  1*1, Fb 


78 
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I  >lu 


DO  528  J=l,6 
OCNl(I,J>=BCON<I, J) 

8CN2< J,I)=BCONT(J,I) 

00.411  NUE 

CONTI NUE 

DO  572  1=1, G 

DO  631  J=i,6 

BCP1 ( I, J) =  BSUP (I ,  J) 

IFCN7.GT.0)  6SPK1 ,  J)  =  8S P*»  (I,  J) 

CONTI  NUE 
CONTI NUE 

I F (NY .  GT ♦ t ) GO  TO  563 
DO  63 k  1*1, F5 
00  533  J=i,N6 
CCM1 (I, J) =CCONT (I , J) 

CCN2< J,I)=CCON(J,I) 

CONTI NUE 

CONTI NUE 

00  536  1=1,  N6 

DO  635  J=1,G 

CSP1( I,J)*CSUP<I, J) 

CONTI NUE 
CONTINUE 

CALL  NRIC(FE,ACN2,SS,Q,P,CX,T0L»I£^) 

Li=N6 
M0=N6 
N0  =  F5 
I  A  =  8 
1 8=8 
IC=8 

CALL  VMULFF (RR1 , CCN2 , LI , MO, N 0 ,1 A ,  I  3  ,  RX , 1 0 , IER) 
M0=F5 
18=16 

CALL  VMULFF  <RX,P, LI, MO,NO,IA ,13, 0<TRN,IC, I EP) 
DO  553  1*1,  F5 
DO  54  C  J=lf  N6 
OK(I, J)=-OKTRN<J,I) 

CONTI NUE 
CONTI NUE 

IF (N7  »GT • 0) GO  TO  671 
DO  67fl  1*1, F5 
DO  660  J=i,F5 
FO(I,  J )  *  £>  •  0 
IF(I, EO.J) FO(I,  J)  * 1«  u 
CONTI NUE 
CONTI NUE 
L1=F5 
M0=6 
N0  =  5 
I A  =16 
18=6 
I C =1 6 

CALL  VMULFF(8C0N,fiRR,Ll,M0,N0,IA,13,NV,i:,IEM 
N0  =  F5 

CALL  VMULFF <RV,BCONT, Li »  MO,  NO » I A  ,  I  3 ,  SO ,  I C ,  IER) 
CALL  MRlC(F5,ACNi,SO,FO,S,ULX,TOl!,IE*l 

Ll=6 


PROGRAM  TETRA 


7k/7k  0PT*1 


rT  N  ‘♦.6*518 


5 15 


520 


5  25 


530 


5  35 


54C 


545 


55r 


555 


5bU 


565 


570 


MO  =  5 
I A  =6 
19  =  5 
I  C=5 
NO  =  c5 

CALL  VMULFF(R*,BCN2,L1,M0,N3 ,IA,  13,  R7  ,  IC  ,  IER  > 

H0*F5 

I9  =  1C 

CALL  VMULFF(KY,S,L1,M0,N0,H  ,IB,  50,  I C ,  I 
571  L1=F5 

M0  =  5 
NO  =  Ff> 

IA  =  16 
19=6 
IC  =  16 

CALL  \)MUL  FF  (BCN1,  GO,  LI,  MO»NO,  I  A,  1 3,  US  ,  IC ,  I ER  > 
M0=N6 
19  =  9 
IA  =  15 

CALL  VMULFF <0K,CCN2,L1, MO, NO, IA,I 3,  C3,1C, ILR) 

00  592  1*1,  F5 
00  691  J=1,F6 

cctci  ,.j)*cc(j,i> 

591  CONTINUE 

592  CONTINUE 

00  69  A  1=FF,F5 

00  593  J=1,N5 

CN (I , J> "CCT (I-N5, J+N5) 

597  CONTINUE 

694  CONTINUE 

00  596  I=FF,Fb 
DO  595  J=FF,F5 
CN(I« J)*CCTCI-N5, J-N5) 

595  CONTINUE 

695  CONTINUE 

DO  698  1*1,  N5 
00  597  J*1,F5 
CN  (I  ,  J)  *<;  .J 

597  CON'.  I NUE 

598  CONTINUE 

00  730  1*1, F5 
DO  728  J=1,F5 
C"'T. J)*CN(I,J) 

CO (X » J)*ACN1 (I,J> -CN (I, J) 

ULd,  JI*ACN1<I,JI  *UC  (I,  J) 

COC (T  ,J)=ULX(I,J) -CX <1,  J)  *ACN2(I,  J> 

720  CONTINUE 
73  C  CONTINUE 
L1*F5 
M0*N6 
N0*5 
IA*16 
19*5 
IC*15 

CALL  VNUIFF(GK,CSP1,L1,HO,NO,IA,I3,03PIIL,XC,IER> 
L1*G  80 


P*0G3AM  TET^A 


74/74 


OPT  =  1 


c  T  N  *♦.  8*i>18 


MO=6 
NO  =  F5 
I  A  =  lb 

575  IC=t6 

CALL  VM’JLFF(BSP1,  GO,  LI,  MO, NO  >IA  ,  1 3 , 0 S PIL  L ,  1C  , IL*> 


0  CREATE  TOTAL  SYSTEM  MATRIX 

58C  0 

no  75  0  1=1, F5 
DO  7 4 0  J=1,F5 
79  ( I ,  J)=ULX(I,J) 

585  7  40  CONTINUE 

75  0  CONTINUE 

DO  *7  0  1=1, F5 
DO  76  0  J  =  Fl,F2 
70 (I,  J)=UC(I , J-F5) 

59C  7  61  CONTINUE 

771!  CONTINUE 

DO  79  1  1=1, F5 
DO  78  9  J=F3,F4 
79(1,  J)=o.L 

595  73  P  CONTINUE 

79  P  CONTINUE 

DO  310  I=Fl , F2 
00  30  0  J=l, FS 
79d,  J!  =0  #  u 

63f  310  CONTINUE 

919  CONTINUE 

DO  83  0  I=F1,F2 
DO  820  J=F1, F2 
79(1,  J)=CX(I-F5, J-F5) 
oC5  321  CONTINUE 

330  CONTINUE 

DO  85  0  I=F1, F2 
00  840  J=F3,F4 
71  C,  J)  =  0SPILL  (I-F5,  J-F2) 
610  349  CONTINUE 


ul? 


62; 


625 


950  CONTINUE 

DO  370  I  =  F3,F4 

DO  860  J=1 »  F5 

79(1, J)*CSPILL(I-F2,  J) 

850  CONTINUE 
87  P  CONTINUE 

00  890  I* F3, F4 
00  880  J=F1,F2 
79(1, J)=CSP1LL(I-F2, J-F5) 
9  9  9  CONTINUE 
990  CONTINUE  • 

DO  910  I*F3,F4 
00  30  0  J«F3,FL 
73(1,  J)*ASUPCI-F2, J-F2) 
909  CONTINUE 
919  CONTINUE 

00  312  1=1, 2t 

81 


•MrflMffiUUilMMai 
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DO  911  J=l,2 it 
73 (I, J)=Z9(1,J> 
CONTI  NIJE 
CONTI NUE 
I JOB= 0 
NN  =  ?"-  E5+G 
14  =  20 
I*  =  2P 


EIGENVALUE  ANALYSIS 


CALL  EIGRF(79,NN,IA,IJ0B,W,7C,I7,,K,I  ER> 

DO  914  1=1,2* 

DO  91 3  J=l,20 
79 (I,  J) =Z8 ( I , J) 

CONTI NUE 
CONTI NUE 
PRINT*,"  " 

PRINT *,"0=",BB 
PRINT*,”  ” 

PRINT*, "THE  SYSTEM  EIGENVALUES” 

PRINT*,"  “ 

PRINT*,”  ” 

PRINT*,”  *,  ”t  IGENVALUE* 

PRINT*,"  - ■ 

NN=NN -1 

DO  93  0  1  =  1, NN, 2 
PRINT*,"  " 

PRINT*,”  ",*<!) 

J=I*1 

PRINT*," 

CONTINUE 
I A  =  1F> 

17=8 

CALL  EIGRF<ULX,F5 ,IA , I JOB , WN , 70D , I  7,  M K, I i R) 
PRINT*,"  “ 

PRINT*,”  “ 

PRINT*, "THE  EIGENVALUES" 

PRINT*, "OF  A  +  B3",  "  -."EIGENVALUE* 

PRINT*,"  - •• 

F5=F5 -1 

00  931  I=1,F5,2 
PRINT*,"  " 

PRINT*,"  ” ,  H  H  (  I ) 

J=I*i 

PRINT*,"  ”  ,  W  H  (  J ) 

CONTI NUE 
F5=F5 «1 
DO  933  1=1, F5 
00  932  J* 1, FB 
CY (I, J)=CX(I, J) 

CONTINUE 

CONTINUE 

CALL  fIGRF (CX, F5 , IA, I JOB, Hi,  71, 17 , HK1 , 1ER) 

DO  93  5  1=1, FB 


' ,  W  ( I ) 


82 


r~ 


PROGRAM  TETRA  74/7>*  OPT=l 


FTN 


685 


69? 


695 


700 


7  0t 


7  if! 


715 


7  20 


725 


7  30 


736 


7<»0 


934 

935 


95  2 


95  5 


967 
96  8 


969 

970 


956 


DO  934  J=1,F5 
CXII, J)=CY(I,J) 

CONTI NUE 
CONTI NUE 
PRINT*,”  ” 

PRINT*,"  ” 

°RINT *, ”THE  EIGENVALUES” 
PRINT*, ”0F  A-KC”,“ 
PRINT*,” 

F5=F5 -1 

DO  952  1=1, F5, 2 
PRINT*,”  ” 

PRINT*,” 

J=I*i 

PRINT*,” 

CONTINUE 
F5=F5  +1 

CALL  FIGRF(C0C,F5 ,IA,IJ09 
PRINT*,”  ” 

PRINT*,”  ” 

PRINT*, "THE  EIGENVALUES” 
PRINT*, "OF  A+BG-KC",” 
PRINT*," 

F5=F5 -1 

00  955  1=1, Fb, 2 
PRINT*,”  ” 

PRINT*,” 

J=I*1 

PRINT*,” 

CONTI  NUE 


”,  ”£  IGE  N  V  AL  UF” 

9 

",K1< J) 

M2, 72,1 7, WK2,IER> 

”, “EIGENVALUE" 

M  M  —  —  — __  *• 

9 

" ,  W2 ( II 
”,W2< J) 


F5  =  F5*1 
00  96  8  1*1,  Ft* 

DO  967  J=1»F5 
ACN6(I,J)=ACN1(I, J) 

CONTI NUE 
CONTINUE 

CALL  EIGRF(ACNl,F5»IA,IJOB,M3»73,I7»WK3,IEk) 
00  970  1  =  2,  F5 
DO  969  J=l,Fb 
ACN:'.  ( I,  J)  =ACN6  (I  ,  J) 

CONTI NUE 
CONTINUE 
PRINT*,”  " 

PRINT*,”  ” 

PRINT*, "THE  EIGENVALUES" 

PRINT*, ”0F  A  EIGENVA.UE" 

PRINT*,”  - “ 

F5=F5  -1 

00  956  I  =  i,F5,2 
PRINT*,”  " 

PRINT*,”  ",M3(I> 

J*I*1 

PRINT*,”  ”,W3(J) 

CONTINUE 
F5*P5 *1 


•  8*-5l8 


83 


P^OG^AM  TETRA 


74/74  OPT  =  1 


FTN  4.  8*518 


3  INITI AL  CONDITIONS 


7  45 


751 


755 


76.; 


765 


Y <  1>  =  -.Oil 
Y (2)  =  .006 
Y<3>= -.003 

Y  (4)  =  .01 

Y  (  5  )  =  ••iil  l 

Y  ( 6)  =  .0u6 

„  Y(7)  =  -.003 
Y(3>=  .01 
Y(9>= .001 

Y  ( 13 )  =-.009 

Y  ( 11 ) = • U  G  6 

Y  ( 12) =-.001 

Y  ( 13 ) =-.002 

Y  (14) =. C u  2 

Y  (15)  =.03 

Y  (15)  =-.02 
Y(  17)  =• 0  2 

Y  ( 13 )  =-.C2 

Y  ( 19) =-.003 

Y  (  20 )  =.004 


INTEGRATE  STATE  EQUATIONS 


77C 


775 


760 


78F 


791 


TI=0. 0 
TOUT  =  .1 
IFLAG  =1 
NEQN= 20 

ARSERR=l.GE-03 
RELERR=i. UE-C3 
J  J=1 

1102  CALL  OOE<XDOT,NEQN,Y,TI,TQUr,kELLU,4  8SERR,IFLAG,WORK,  IWORO 
Ic ( TI  .LT .TOUT) GO  TO  11S2 
XX1( J J)=Y (1) 

XX2(J J)=Y(2) 

XX3UJ)=Y(9) 

XX4(  J  J)  =Y  ( If ) 

XX5( J  J)=Y (11) 

XX5C  J  J)*Y  (12) 

XX7 ( J  J)=Y (13) 

XX8(  J  J)=Y  (ID 

J  J=JJ  41 

TOUT=TOUT4.1 

IF (TI «LE. 23 • ) GO  TO  1102 


GENERATE  AND  PLOT  LZNE  OF  SIGHT  X  AND  Y 


795 


DO  1103  1=1. 2C1 
THCI)  «(I-l)*.l 
119*  CONTI  HUE 

00  1104  I*l»  201 

84 


imm mmSBu 


? 900? AM  TETra 


74/71*  OPT  =  1 


FTN  4.8*518 


8Cf’ 


605 


61.' 


8  ID 


8  2.: 


62i. 


8  3C 


835 


841 


845 


85*1 


119 

110) 


11?’ 


957 


958 

959 


96  r 

961 

962 


963 

954 


965 

9S6 


XI (I) =-,2470727*XXl(I)+. 3998 955* XX2CI) 

XI (I) =X1(I)-.087  8  33J1+XX3(I> +.0535/94 4+ XX 4 (I) 

XI (I >  =X1(I) +  .C274  558  6+XX5CI) -.026614C1*  XXL  (1 ) 

XI  (I)  =Xl(I)-.C299  367t  +  XX7CI)  +  . 0 0U& 1 3 5322 3 *  XX 8 < I ) 

X2(I) =.4278565+ XX 1(1) ♦. 2 3 0 529 1+XX 2 ( I ) 

X2  ( I )  =X2<I)-.<i507  0142+XX3<I)  +  .i*35/r784+XX4<I) 

X2  (I)  =X2(I) -. 0475702 2+XX5 CI>  +  .^4bD'iS5  3*'XX6(I) 
X2(I>=X2(I)-.Gi7  3!J931+XX7(I>  +  .  435 1  i  ui.  CD  0 1 219  161  *XX8  (  T) 
CONTI NU£ 

CALL  PLOTU.,0.,-3> 

CALL  SCALE(TM,8.  , 201,1) 

CALL  SCALE(X1,8., 201,1) 

CALL  AXIS(U.,C.,4HTIME,-4,8.  ,0.  ,  1 1  (  24  2)  ,  V  M  C2C3)  ) 

CALL  AXIS  (ti • ,C  «  ,  4HL0SX,  4 , 8.,90.»X1(232),<1  (2  C3) ) 

CALL  LINECTK, XI, 201,1,5,2) 

CALL  SYMBOL  (4.,6.  ,  0. 21,  13HL0SX  VS.  TIME,  C.  ,13) 

CALL  PLOT  (If;.,  0., -3) 

CALL  SCALE(TM,6. ,201,1) 

CALL  SCALE(X2,8. , 201,1) 

CALL  AXISCU.,0.,4HTIME,-4,8.  ,  d.,TK2u  2)  ,  TM  (2!! 3)  ) 

CALL  AXISCC.,C.,4HL0SY,4,8.,9u.,X2(2G2)  ,  <2(203) ) 

CALL  LIN£<TM,X2,20l,l»5»2) 

CALL  SYMBOL  (4., 6., 0.21,  13HLDSY  VS.  TIME, 0. ,13) 

CALL  PLOT E ( N ) 

IF ( M. LT *N7>  GO  TO  1101 

I F  ( <K  ,GT . w ) GO  TO  1101 

05=2. 0 

Q5=.0  01 

00  959  1  =  1, M 

DO  95  8  J=1,12 

PHIT1  (I, J)  =  b.O 

CONTI NUE 

CONTINUE 

DO  96  1  I  =  GG,G 

DO  960  J=l,  12 

PMIT1  (Z,J)=PHIT  (I  +F5  ,  J) 

CONTINUE 
CONTINUE 
L1  =  G 
MD  =  12 
no=3 
I A =1 6 
19=12 
I  C  =  16 

CALL  VMULFF(PHIT1,D,L1,MO,NO,IA,IB,9SP3,IC,IER> 

00  964  1=1, G 
00  963  J=  1,  t> 

BSP4C I,J)=3SP3(I,J) 

CONTI NUE 
CONTINUE 
DO  966  1  =  1, G 
00  96  5  J=1,G 
C(I,J)«0.t 
I*CI« EO. J) C (I , J) =1 .0 
CONTINUE 
CONTINUE 
IA  =  16 


855 


85 


861 


PROGRAM  TFTRA  7A/74 


M0  =  3 
N0=6 
IC  =15 
N9=G 


OPT  =  1 


FTM 


865 


871 


!  875 


CHECK  FOR  ZERO  SINGULAR  VALJL 


CALL  LSVDF(6SP3,IA,N0,N0,C,IC,N8,S,  WK,IER> 
Q4=S ( 6> 

1^(94. LT.Q5)GO  TO  1075 
QS  =08 -.81 

IF  (16  .LT.-2.OGO  TO  1C96 
01  =  ( 4  .3+Q6**  2,)** .6 
D(7,3)=l./01 
0(8,3  >=1.7318/01 
0  0,3  >  =Q6/Q1 
GO  TO  962 
1)7'  PRINT*,"  " 

PRINT*, "THE  LEAST  SINGULAR  i/ALUE  IS  ",S(5> 
PRINT *,"06  =  ", Q6 


880 


APPLY  TRANSFORMATION  TECHNI3UE 


885 


89C 


895 


9  JO 


I  9')6 

I 


I 

i. 


[  9 1C 


00  1076  J  =  l>6 
TTC1,  J)  =  6SP3 (C, J) 

1176  CONTINUE 

no  10 78  1=1,6 
T ( I , 1 ) =TT (1,1) 

1978  CONTINUE 

00  1889  1=1, N5 
00  10  79  J=l,12 
PHIT2  (I,J)=(/.0 
1970  CONTINUE 
118  CONTINUE 

DO  10  82  I=FF,F? 

00  1081  J=l, 12 

PH  IT  2  (I,J)=PHIT (I ♦M,  J) 

1981  CONTINUE 
198?  CONTINUE 
L1=F5 
M0  =  12 
N0  =  6 
I  A  =  16 
19  =  12 
IC  =  16 

CALL  VMULFF (PM1T2 ,0,  LI,  MO  ,  NO ,IA,  19,  SCON, 13 ,1ER) 
M0=5 
N0  =  1 
19  =  5 

CALL  VMULFF ( SCON, T, LI, MO, NO, IA , 1 9 , 9 T , IC, I E R) 

L  1*1 
N0*6 
I  A»1 


. 8+518 


I 

! 

i 

I 

I 

! 

i 


86 


PROGRAM  TFT  R  A 


7A/74  OPT  =  1 


FT  N  4. 8+518 


y 

915 


9  20 


925 


93‘3 


935 


94'j 


945 


950 


955 


96( 


965 


IC  =  1 

CALL  VMULFF(TT,R,L1,M0,N0,IA ,18,  RTl  ,  IC  ,  I E  R  > 

NO  =  1 

CALL  VMULFF(RT1,T ,ll,M0,N0,IA ,13,Rr,lC,lER) 

RTI (1 ,1)  =  1./RT (1, 1) 

L1=F5 
40=1 
IA  =  lo 
19  =  1 
i  c = 1 6 

CALL  VMULFF(BT,RTI,L1,M0,N0,IA,I3,  VI,  1C,  IER) 

DO  l^S  J  =  1,F5 
9TT(1,J)=BT(J,1) 

13©-  CONTTNUE 
L1  =  F5 
M0  =  1 
N0  =  F5 

r  a =16 
18  =  1 
IC  =  lb 

CALL  VMULFF (V1,BTT , LI, MO, NO, IA, IB, 77, IC, IER) 

DO  i'J  86  1=1, F5 
00  10  85  J  =  1,F5 
QO { I , J)=u .0 

IF (I. EO.J)QC<I, J) =1000.0 
U8r:  CONTINUE 

10  8r  CONTINUE 

CALL  MRIC (Ft  » ACN1 ,  VV  ,  QQ  ,  SI, JLX1, T  DL , ItR) 

L  1  =  1 
40  =  1 
NI0  =  F5 
I A  =1 
19  =  1 
IC=1 

CALL  VMULFF (RTI,BTT,L1,MO,NO,IA,I3,  301 ,13, IEP) 
M0=F5 

n=tb 

CALL  VMULFF  < G01 , SI, LI, HO, NO,  IA, 13, 302 , 1C, IER) 

L 1 =6 
M0=1 

N0="5  ‘ 

I A =6 
18=1 
I  C=5 

CALL  VMULFF (T, GO 2, LI, MO, NO, IA, IB, 53, IC, IER) 

00  1086  1=1,6 
00  13 87  J=l,6 
3SP1<I,JI=BSPMI, J> 

LT8*  CONTINUE 
HA’  CONTINUE 

no  It)  90  1  =  1, F5 
00  ID  89  J*l,6 
9CN1(I,J)=BC0N(I,J) 

1189  CONTINUE 
1 09r  CONTINUE 

00  lit  9?  1  =  1,  F5 
DO  10  91  J*1,F5 

87 


i _ 


AD-A094  766  AIR  FORCE  INST  OF  TECH  WRI6HT-PATTERSON  AFB  OH  SCHOO— ETC  F/6  12/2 

MODERN  OPTIMAL  CONTROL  METHODS  APPLlEO  IN  ACTIVE  CONTROL  OF  A  T-- ETCtU) 
DEC  00  A  M  JANISZEWSKI 

UNCLASSIFIED  AFIT/GS/AA/0OD-2  NL 


PROGRAM  TETRA  74/74  0PT  =  1  FTN  -4. 0*516 


971  ULX(I , J)=ULX1CI, J) 

1191  CONTINUE 
L  J  9  ’  CONTINUE 
KK=KK <1 
DO  1201)  1=7,9 
975  D1 <5, I)=D(I,3) 

120  CONTINUE 
N7  =  l 

GO  TO  118 

1 19f  PRINT*, -SINGULARITY  PROGRAM  FAILE  D" 
981’  1101  STOP 

END 


If 


It 


21 


Zb 


3f 


35 


40 


45 


55 


55 


SUBROUTINE  MRIC(N,A,S,Q, X,7,TOL,IER) 
DIMENSION  A<16) ,S<26) ,Q<16) ,X (15) ,7(16) 
C0MMUN/MAIN1/NDIM,ND1M1 , F(1S) 
C0MM0N/MAIN2/TR<16) 

COHMON/INOUT/NOT 
AOV=TOL*1  .E-t'6 
NN  =  N*  NOIM 
NM1=N -1 

I  NO=l 
COUNT =0 • 

IE (IER.EQ.l) COUNT =99 . 

IF(IER.EQ.1)MR=N 
IF  (IER.EG.i)  60  TO  10  u 
ri=-i . 

300  CONTINUE 
IER=U 

COUNT  =COUNT+l. 

DO  15  1=1, N 
00  15  J  =  I  ,  NN,  NDI  M 
15  X(J)  =  -S(J) 

CALL  INTEG( N»A, X  »  Z,T 1) 

CALL  FACT0R1(N,Z,X,MR) 

IER=i 

IE ( MR  *LT . U ) GO  TO  20 0 
IER=0 

CALL  GMINV(N,N,X,Z,MR,0) 

CALL  TER(TR,Z,N,N,1,  2) 

CALL  MMUL ( Z,TR,N , N ,N , X) 

DO  18  11=1, NN,NOI Ml 
1  =  11 

DO  17  J=II,NN,NOIM 
X  (J)=  (X(J)+X  (I)  >  /  2% 

X(I)  =  X(J) 

17  1=1*1 

18  CONTINUE 

101  CONTINUE 

DO  16  1=1, N 
15  TR(I>=-1.0 

T0L1=  TOL/lli  • 

MAXIT  *4C 

DO  40  IT=l,hAXIT 
IF (lER.EQ.l) GO  TO  10 i 
CALL  MMUL (S,X,N,N,N,  F) 

CALL  MMUL (X,F,N,N,N,  Z) 

DO  20  1=1, NN, NOIM 

II  =  1  +  NM1 

DO  20  J«I,II 
X( J)= A(J) -F( J) 

20  7(J)  =  7(J)  *0  (  J ) 

101  CONTINUE 
IER*0 

CALL  MLINEQ(N,X,7,X,T0L1,IER) 

I E  ( I E  R,NE  •  U )  GO  TO  200 
L  =  0 

C1=0. 0 

11=1 

00  25  1=1, N 


{ 

!  i 

■A 

■1 


89 


SUBROUTINE  MRiC 


74/74  OPT  =  1 


FTN  4.8*518 


bG 


bb 


7C 


7b 


8? 


85 


IF(ARS(X(II)-TR(I>>.LT. ( AOV*TOL »X ( 1 1 )  ) >1  =  1*1 
TR  ( I ) =X(II> 

II  =  II  ♦NOIM1 
25  Ci=Cl *TP ( I ) 

IF(ABS(Ci)  .GT.1.E*2C)G0  TO  5  o 

I F (L • NE.N) GO  TO  40 

CALL  GMINV(N,N,Z,F,MR,0) 

CALL  MMUL (S,X,N,N,N, Z) 

DO  3C  I=l,NN;NOIM 
II  =  I*  NM1 
DO  30  J=I,II 
33  Z( J)=A(J)-Z(J) 

I F  <  MR  .NE.N) WRITE (NOT ,35) MR 
35  FORMA  T(26H0RICCATI  SOLN  IS  »SD--RANKl 3 ) 

GO  TO  65 
40  CONTINUE 

WRITE  (NOT, 45) MAXIT 

45  F0RMAT(2fcHfcRICCATI  NON-C ONVERGENT  INi2,llH  ITERATIONS) 

GO  TO  63 

5.)  WRITE  (NOT, 55)IT,T1 

55  FORMAT ( 29HGR1CCAT I  BLOW  UP  AT  ITERATI ONI  2, 12H  INITIAL  T=F13.5) 

63  IER=1 

65  RETURN 

23«  IF  (IND.EQ.2)  GO  TO  250 
I F (COUNT. GE. 1C.) RETURN 
T1  =  T1 / ( 2, **C  OUNT ) 

IN  0=  2 
GO  TO  300 

250  Ti=Ti*(2.**C0UNT> 

I  NO= 1 
ENO 


FUNCTION  XNORM  74/74  OPT  =  l 


FTN  4.8*518 


1 


5 


10 


15 


2r 


13 


15 

23 


FUNCTION  XNORM (N, A) 
DIMENSION  A (16) 

COMMO  N/MAIN1/N0IM, NO I Ml 

NN  =  M*  NOIM 

01=0. 

TR=A (1) 

IFtM.EQ.DGO  TO  20 
1  =  2 

DO  10  II=NOIN1,NN,NOIM 
J  =  II 

00  5  JJ=I ,11, NOIM 
C1=C1 *ABS (A ( J)*A ( JJ) ) 

J=  J*1 

TR=TR  *A( J) 

1  =  1*1 


TR=TR /FLOAT (N) 

DO  15  II*1,NN,N0IM1 
C1=C1 *(A(II)-TR)**2 
XNORM  =ABS(TR>  *SQRT (Cl) 
RETURN 
ENO 


# 


90 


n 


SUBROUT  INF  *1LT  NEQ 


T  4  /  7  4  OPT  =  i 


FTN  4.8*518 


SUBROUTINE  MLINEQ(N,A,C,  X,TOL,IERI 
DIMENSION  A  < 16) , C ( 16 )  »X  ( 16) 

C0MM0N/MAIN1/NDIM,ND1M1 
COMMON/MAIN3/F (16) 

AOV=T  OL*i.E-06 
DT  =  .  5 
OTl=U  . 

NN=N*  SOI M 

HO  5  II=i,NN,NOIMl 
DT1=DT1-A (13 ) 

OT 1  =  0  71/N 

IF<OT1.GT.4.1.>DT=OT*4.P/OT1 
TI  =  1 

DO  2C  1=1, N 

DO  15  JJ=I,NN,NOIM 

X( JJ) =OT»A(JJ> 

X(II) =X(II) -.6 
II  =  II <N0IM1 

CALL  GMINV(N,N,X,F,MR,U) 

IER=4 

IF(MR  .NE.N) RETURN 
CALL  MMUL <C,F,N, N,N, X) 

1  =  1 

DO  46  11  =  1 , NN» NDI M 
J=II 

IF(I.  EO.  1) GO  TO  3 J 
DO  25  JJ= I , II , NDI M 
C(J)  =  C(JJ) 

J=J* 1 
1 0  =J 

00  35  JJ=II  , NN,NOIfi 
C ( J)  =  CT*DOT (N,F(2I),X(JJ>) 

J= J*1 

F  ( ID)  =F(ID)+l.t 
1  =  1*1 

DO  90  IT=1,20 
NEZ=0 

CALL  MMUL (C,F,N,N,N,X) 

1  =  1 
11=1 
J=1 

GO  TO  70 
J*I’ 

DO  65  JJ=1, II, NDI H 
C< J)  =  C(JJ) 

J=J*i 
10  =  J 

OT 1*C  (J) 

no  75  J J= II , NN, NOIM 

C(  J)  =  C(J)  ♦OOTCN,  F(II)  ,X(  JJ) ) 

J=  J*i 
J=J-1 

DO  80  JJ= II ,  J 
X( JJ>  =F ( JJ) 

IF(ARS<C<ID)  ).GT.1.£*15J)G0  TO  95 

IF  MB  SCC  (10)  -DTI  I  .LT.(ADV*TDL*A3S(I(ID) ) )  )  NE7=N£7*1 
1*1*1 


5 UBRD UT I NF  MLTNEO 


74/74  OPT  =  1 


FTN  4.8+518 


II=II +  NDIM 

IF (I , L£«N) GO  TO  60 

F*».i 

IF(NE,.ED.N)GO  to 

15  C 

CALL  MMUL (X,X,N,  N, 

N,  F) 

TO 

CONTI  NIJE 

95 

IER=i 

PET UR  N 

6f 

15  ? 

CONTI NUE 

NM1=N -1 

DO  155  1=1, NN, NO  I M 
1 1  =  I  *  NH.i 

DO  155  JJ=I,II 

7i 

155 

X ( JJ) =C( JJ) 

IER=C 

RETURN 

END 

SUPRDUTTNF  FA^TOPl  7*./74  CPT=i 


FTN  <♦.  8+518 


1 


r 


1- 


15 


2u 


25 


si- 


subroutine  FACTORl(N,A,S,MR) 

DIMENSION  A (It ) »  S ( 16 ) 

COMMON/MAIN1/NDIM,NOIM1 
COMMON/INOUT/KOUT 
TOL  =  l  .E-l>6 
MR  =  0 

NN=N* NDIM 
TOLi=P. 

DO  1  I=1,NN,NDIM1 
R  =  A9S  (A  (I ) ) 

l  IF(R.GT.TOLl)TOLl=R 

TOLl=TOLl*l.E-12 
11  =  1 

DO  50  1=1, N 
I M 1=1  -1 

DO  5  JJ=I  ,  '.»N,  NDI M 
5  S( JJ)  =0. 

10  =  11  +IM1. 

R=A (I D) -DOT (IM1,S(II),S(II)> 

IF(A5S(R)  .LT.  (TOL*AUD> +TOL1  ))GO  ID  30 
IF(R> 15,51, 2C 
15  MR  =  -i 

WRITE  (KOUT, 1 juC ) 

im  F  0  RFA  T  ( 37  HGTR.l  ED  TO  f  ACT  Ofl  AN  INDEFINITE  MAT  FIX) 

RET JR  N 

20  S (ID) =SQRT (R) 

MR=MR+i 

I F ( I , EO.N) RETURN 
L=  II*  KOI M 

00  25  JJ=L, NN, NOI M 
I  J=JJ*IM1 

25  S(IJ) ■(A(ZJ)-OOT (IH1,S(ZI),S( JJ) ))/S(IO) 

5)  IT=II*NOIM 

return 
END 


35 


92 


SUBROUTINE'  INTEG 


74/74  OPT=i 


FT N  4. ft+518 


10 


15 


Zl 


2‘i 


3u 


35 


4C 


45 


6C 


55 


SUBROUTINE  I  NT  EG  (  N  ,A  , C»  S  ,  T  ) 
DIMENSION  A<  16)  »  C  (16)  »S  ( 16) 

C0MM0N/MA1N1/NDIM,N0IM1 , X(15) 

COMMO  N/MAIN2/C0EF (16 ) 

NN  =  N*  NDIM 
NM1=N -1 
I  ND= ') 

ANORM=XNORM (N*A) 

DT  =  T 

IF(AN0RM*A9S(DT)  .LE.fc.5)  GO  TO  10 
DT  =  OT /2 • 

IND=I NO  +  1 
GO  TO  5 

1!)  DO  15  1=1, NN, NDIM 

J  =  I«-NMl 
DO  15  JJ=I,J 
15  5 ( JJ) =OT*  C( JJ) 

Tl=0T**2/2. 

DO  25  IT=3, 15 

CALL  KMUL (A ,C,N, N ,N, X) 

03  20  1  =  1, N 
II=(I -l)*NOIM 
DO  2C  JJ=I , NN, NDI M 
II=II *1 

C(JJ)=(X(JJ)+X(II))*T1 
20  S( JJ)  =S(JJ) *C(JJ) 

25  T1=DT /FLOAT (IT) 

IF(IN0.E0.o)G0  TO  1J 0 
COEF(  11)=1.(» 

DO  30  1  =  1,11’ 

11=11-1 

TO  CO£F( II) =DT*COEF ( II *1) /FLOAT (1) 

11=1 

03  40  1=1, NN, NDIM 

J=I*NM1 

DO  ?5  JJ=I,J 

35  X(JJ) =A ( JJ) * COEF ( 1) 

X(II) =X (II) ♦COEF (2) 

40  II=II ♦NOIMl 

DO  35  L=3,ll 
CALL  MMUL (A,X,N,N,N,C) 

11=1 

Tl=CO EF(L) 

DO  >5  1=1, NN, NDI M 
J  =  I  *-N  Ml 
DO  50  JJ=I,J 
50  X(JJ)=C(JJ> 

55  II=II  ♦NDIM1 

L  =  0 

50  L=L+1 

CALL  MMUL (X,S,N,N,N, C) 

11*1 

DO  9 ii  1*1, N 
J*II 

IF ( I • EO,l)GO  TO  75 
DO  TU  JJ*I,II,NOIM 
S(  JJ)  *S(J) 

9  3 


SUBROUTING  I  NT  G  G  7L/7  +  OPT=i 


f  r  s  > 


6t' 


b5> 


7 


71  J=  J  +  i 

75  00  95  JJ=I,N 

KK  =  JJ 

DO  90  K=I,NN,NDIM 
S(  J)  =  S(J>  +C(K)+X  (KK> 

91  KK  =  KK  +  NDI M 

15  J=J  +  NOIM 

DO  87  JJ=1 , NN» NOI M 
97  C(JJ)=X(JJ) 

31  I I=II +NDI M 

IF(L. EQ»IND> GO  TO  IOC 
CALL  NMUL (U,C,N,N,N, X) 
50  TO  60 
11"  CONTINUE 
RETURN 
END 


SUBROUTINE  NMUL 


74/7»*  OPT  =  1 


rT  N  4 


1 


c 


10 


SUBROUTINE  KMUL (X  ,Y, N1,N 2, M3 , 7) 

DIMENSION  X  (16)  ,Y  (16)  ,Z  (16) 

COMNO  N/MAIN1/NDIM 

NGN03=NDIM*N3 

NEND2  =NDIM+N2 

DO  1  1=1, Ni 

00  1  J=I»NEN03»NDIM 

7( J)=  9,0 

<<  =  J-  I 

DO  1  K=I, NEN02,N0IM 

7*}  J..X  ««.,<«» 

RETURN 

ENO 


FUNCTION  f)0T  74774  OPT=l  5TN 


1  FUNCTION  OOT  (NR, A , B) 

DIMENSION  A  ( 1 6 )  ,  B  ( If  * 
f)OV=0  . 

DO  1  1=1, NR 

'  1  DOT  =  COT+A (I )*0 ( I ) 

RETURN 

ENO 


, 8+518 


, 8+518 


4. P+518 


94 


SUBROUTINE  T  Ft 


74/74  0PT  =  1 


FT  N  4. 8*518 


1 


1C 


15 


20 


2'j 


3f- 


35 


40 


S'JDROITINE  TFR(X,A,N,M,  K,I> 
01  MANSION  X ( 1 6 ) ,  A  (16) 
CD’HON/MAINl/NDIM 
JS=(K-1)*NDIM*M 
JEND=N*NDIM 

90  TO  (10,30,54,7  J, 90), I 
13  DO  20  11=1, N 

00  20  JJ=II,JEND,NDIM 
21  X ( JJ>  =A( JJ*JS> 

RETURN 

31  DO  '♦!)  11  =  1, N 

<<=(I  I-D^NDIM 
DO  40  JJ=1,M 
LL=( J J-l) *NOlM*II 
41  X(<K*JJ>=A<IL+JS) 

RETURN 
51)  KK  =  0 

DO  6>)  11= 1, JENO,NOiM 
LL  =11  +N-1 

no  o*J  J  J=  II ,  LL 
KK=KK *1 

51  X  (  KK)  =  A  ( J  JOS) 

RETURN 

M  KK=M*N*1 

DO  80  11=1, M 
LL=(M-II)*NDIM*1 
DO  50  I J=1 » N 
KK=KK -1 
J  J=LL*N-I  J 
31  A  (  JJ*  JS)  =X  (KK) 

RETURN 

90  SA  VE=  A ( 1) 

K=N 

DO  91  1*1, N 
L  =  NI 

DO  92  J= 1 , N 
IK=(*-l>*NOIM*K 
X ( IK) =0. 

I F (L  « EQ.K)X(1K)=A(L) 

92  L=L-1 

91  K=K-1 

X  ( 1)  =  SAVE 

RETURN 

END 


9  5 


SUBR3UTINF  GMTNV 


7*4/74  0PT  =  1 


PTM  4.8+516 


1 


r> 


1'.' 


15 


2»: 


25 


3L 


35 


43 


45 


50 


55 


SU5R0  UTINE  GMINV  (  NR,  NC,  A , U, MR, MT> 
01  MEN  SION  A(1£),U<16> 

COMMO  N/MA IM1/NPI M , NDI Ml ,  S  (15  > 

COMMON/INOUT/NOT 

T0L=1 • E-l 2 

MR  =  NC 

N®  Ml  =  NR-1 

TOL 1  =  1.  £-2<. 

J  )  =  1 

00  11.0  J=1,NC 
FAC=nOT(NR,A(JJ) ,A(JJ)) 

JM1=J-1 
JRM= J J+NRMi 
JCM=J J+JM1 
00  20  1  =  J J»  JCM 
21  U  (  T)  =  0. 

0 ( JCM  1  =  1. 0 

Ic (J. £0. 1 ) GO  TO  54 

KK  =  1 

DO  3  it  K  =  1 ,  JM1 

IF  (S  (  K)  •  EQ.l  *(.)  GO  TO  30 

T F MP  =  -DOT  (NR,  A  (JJ)  ,A (KK) > 

CALL  VADD(K,TEMP,U(JJ),U(KK)  > 

V)  KK=KK+NDIM 

00  50  L=l,2 
KK  =  1 

DO  53  K=1,JM1 
I ^ ( 3 ( K) . EQ • 6 . ) GO  TO  56 
TEM®= -OOT (NR, A ( JJ) ,A (KK) ) 

CALL  VADD (NR , TEMP , A ( J J) ,A(KK) ) 
CALL  VADD(K»TEMP,U(JJ)  »  U  ( KK)  ) 

5  1  KK=KK  +NDI M 

TOLl=TOL*FAC 
FAC=OOT(NR,A(JJ> , A  (JJ)) 

54  IF(FA  C.GT.TOLDGO  TO  71 
00  55  I=JJ,JRM 

55  A(I)  =  0. 

S(J)  =  3. 

KK  =  i 

00  55  K=i,JMl 
I F ( S ( X) *EQ. 0 • ) GO  TO  65 
TEMP=  -DOT  (K,U(KK)  ,U(  JJ)  ) 

CALL  VADD  (NR, TEMP, A(  JJ)  ,A(«)  ) 

65  k<=KK+NDIM 

FAC’OCT(J»U(JJ) , U ( JJ - ) 

MR  =  MR-1 
GO  TO  75 

71  S ( J ) = 1*6 
KK  =  1 

00  72  K* 1, JMi 
I F  (S  (  K)  .  EQ.  1 )  GO  TO  72 
TEMP* -DOT (NR, AC JJ) ,A (KK) ) 

CALL  VA DOCK, TEMP, U(JJ),U(KK)  ) 

72  KKaKK+NOIM 

75  FAC=1 ./SORT (FAC) 

Of)  80  I*J J, JRM 
81  A(I)*A(I)+FAC  96 


SUBROUTINE  OMfNV  74/74  0PT=1 


PTN 


60 


6r 


70 


7? 


DO  55  I=JJ,JCM 
55  U(I)  =  U<I>*FAC 
10  A  JJ=JJ+NOIM 

IOMR  .EO.NR.OR.MR.EQ.NC)  GO  TO  12S 
T E ( MT  .NE. G)WRITE(NOT,110)NR,NC,MR 
tin  F0RMAT(I3,1HX,I2, 8H  M  RANK, 12) 
12?  NEND=NC*NOIM 
J  J  =  i 

DO  135  J=1,NC 
DO  12*  1=1, NR 
II  =  T- J 
S(I)  =  3. 

DO  125  KK=JJ,NEND,NOIM 
125  S  <  I )  =  S(I)  +A(11+KK)*U(KK) 

1 1  =  J 

DO  133  1=1, NR 
U  (  TI ) =S  ( I  ) 

I3rt  II  =  II+NDIM 

135  J  )  =  J J  +  NDI Ml 
RETURN 
END 


SUBROUTINE  VA  ">  D  74/74  0PT  =  1 


FTN 


1  SUBROUTINE  VADO < N  ,  Cl , A , B ) 

DIMENSION  A ( 16 ) , B  ( 16 ) 

DO  1  1=1, N 

1  A(I)  =  A(I)+C1*B(I) 

5  RETURN 

END 


SUBROUTINE  XOOT  74/74  0PT=1 


FTN 


1  SUBROUTINE  XOOT ( TI  ,Y , YP ) 

DIMENSION  Y  (20  ,YP(20) 

COMMON  Z9<2C,2l>) 

L 1  =  2  3 

5  M0=2D 

N0  =  1 
I A  =  L  ll 
IB  =  20 
IC=23 

1C  CALL  VMULFF (79,Y,L1,M0,N0,TA ,I3,Y1 * 3 * 5,I3,IE3> 

RETUR  N 
END 


4. 8*510 


4. 8+518 


4. 8+518 
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parison  is  obtained  by  complete  elimination  of  the  spillover 
effect.  Using  singi.  ar  value  decomposition,  the  spillover  is 
first  eliminated  through  judicious  reorientation  of  one  sensor/ 
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demonstrates  the  advantages  of  spillover  elimination,  but  fails 
to  satisfy  the  specified  bandwidth  for  error. 

Sensors  are  added  to  the  model  at  the  fourth  mass  and  ob¬ 
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initial  sensor /actuator  pairs  is  no  longer  applied.  Line  of 
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